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1. Introduction. This paper is concerned with the simplest equation of mixed 
type, known as the Tricomi equation 


T(u) = zu,. + u,, = 0 (1) 


where u is the dependent variable, z and z the independent variables. Equation (1) is 
elliptic when z > 0, hyperbolic when z < 0. In this latter half plane the characteristics 
of (1) are the lines defined by 3x + 2(— z)*”* = constant. The Tricomi problem consists 
of solving the equation T(u) = 0 in a domain A bounded by two concurrent character- 
istics AC and BC drawn in z < 0 and an arc AMB drawn in the half plane z > 0 when 
values of u are known along AMB and AC. Thus the mixed character of the equation 
is involved in the definition of the Tricomi problem and in fact the Tricomi problem is 
the typical problem for an equation of mixed type. The existence of a solution for such 
a problem has been proved by F. Tricomi [1]** himself in his fundamental paper. A 
quite different type of proof has been given by P. Germain and R. Bader [2], [3]. They 
have considered in particular the special case for which AMB is a “normal’’ curve, 
according to Tricomi’s terminology—that is to say an arc defined by (x — 2)” + y? = R’, 
3y = 22°”, with x and R given constants and have shown that in such a case it is possible 
to give an explicit solution of the Tricomi problem. Such a Tricomi problem will be 
called a ‘‘normal”’ Tricomi problem. This result was quite interesting from a theoretical 
point of view and permitted a very simple proof of the existence of the solution of the 
Tricomi problem to be given. More recently [4], [5], and [6] it was shown that such a 
“normal” problem was of particular interest in the application of an approximate method 
to subsonic and transonic flows involving jets or wedges. However, the solution given 
previously was not found quite satisfactory from the computational point of view. 
It was emphasized, [2], that a transformation involving three parameters (transformation 
related to the Poincare’s geometry) can be associated with the Tricomi equation. As a 
result, the solution of a “normal” problem can be simply derived from the special case 
in which A is the region A, defined by x > 0 for z > 0, 32 > 2(— z)*” for z < 0. The 
values of u along oz (x > 0) are known; w is also given either along AC, or along the 


*Received May 5, 1955. The results presented in this paper were obtained in the course of research 
sponsored by the National Advisory Committee for Aeronautics under Contract NAW-6323 while the 
author was Visiting Professor of Applied Mathematics at Brown University. 

**Numbers in brackets refer to the bibliography at the end of the paper. 
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characteristic at infinity; in the former case, which will be called the direct problem, 
the solution must be regular at infinity [2]; in the latter case it will be called the conjugate 
problem. An inversion with center at the origin allows one to reduce one of these problems 
to the other. A few comments are needed on the definition of the Green’s function of a 
Tricomi problem. This notion was introduced in [2], [3] for the case of a Tricomi equa- 
tion and generalized in [7], [8] for a wider class of equations of mixed type. It arises when 
one looks, for the possibility of writing the expression for the solution of a Tricomi 
problem as a linear functional of the data. It can be shown that to every point P inside 
A, one can make correspond a function gp(M), called the Green’s function of the Tricomi 
problem for the given domain A. This function is continuous for M inside A, (if some 
singular lines or points are excluded), and has the following fundamental properties: 
Given one solution u of (1) in A, it is possible to find by application of the Green’s 
formula the value of u(P) in terms of some integrals involving only (besides the value’ 
of gp and its derivatives), the value of u along AMB and AC. Moreover, gp as function 
of the coordinate of M is a fundamental solution of (1) and gp is zero on AMB and BC; 
gp as a function of the coordinates of P, M being kept fixed inside A, is a fundamental 
solution of (1) and takes the value zero when P is on AMB or on AC. In other words, 
as a function of M, gp(M) satisfies some boundary conditions for the conjugate problem; 
as a function of P, it satisfies some boundary conditions for the direct problem. The 
notion of fundamental solution* for an equation of mixed type is also discussed in [7], 
[8]. When the point P is in the elliptic half plane (z > 0), gp as a function of M is regular 
everywhere in the open domain A, except in the neighborhood of P, near which it has 
the classical logarithmic singularity. When the point P is in the hyperbolic half plane 








Fia. 1. 


(2 < 0), the Green’s function as function of M is singular on the characteristics shown 
in Fig. 1: it has some discontinuities along PQ and PS, proportional to the values of the 
Riemann function along these lines, and it becomes logarithmically infinite in the 
neighborhood of the reflected characteristic QR. A similar behavior is of course valid 
for g as a function of P when M is fixed inside A. 


*In the usual French terminology such a solution is called “elementary”. 
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The following developments are the result of two remarks. First it was shown in 
[7] and [8] how it is possible to build the Green’s function of a strip for a class of differ- 
ential equations, even if the strip lies in a mixed region. Second, for the Tricomi equation 
(1) new independent variables can be introduced in such a way that the Green’s function 
of the Tricomi problem corresponding to Fig. 2 is transformed into the Green’s function 
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of a strip in a mixed domain. Thus it is possible to apply a previous technique with 


some minor differences in order to obtain the required result. 
2. Transformations of the fundamental equation. The following definitions were 


introduced in [2]: 
By= 2°", rP=ert+y, r=rt. (2) 
In A, , r’ as defined by (2) is positive (zero on the characteristic OC), and accordingly 
r will be assumed real and positive. With r and ¢ as new independent variables, the 
operator T'(u) becomes: 
T(u) = 2L(u) = z{u,, + (1 — &)r*un, + (4/3)¢"'u, — r*tu,)}. (3) 
Now it is convenient to introduce the new variable ¢ defined by r = exp ¢, and with 
€ and ¢ as independent variables (3) becomes 
T(u) = zL(u) = r-*zM(u) = r7*z{(1 — Our + ge + 1/8(u, — 4tu,)}. (4) 
The domain A, is mapped into the half plane ¢ > 0 of the é, ¢ plane. For 0 < ¢ < 1, the 
equation M(u) = 0 is of the elliptic type, and for ¢ > 1, hyperbolic. Another form can 
be used with £ and ) as independent variables where 


r= fia-re, 


namely 
T(u) = r-2M(u) = r72(1 — #7)? N(u) 


r2(1 — #)-? fun + (1 — #)'(uge + (1/8)up} ©) 
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The Green’s function we are looking for is a fundamental solution of (1); precisely, such 
a fundamental solution ep is a solution of 


T(u) = 8... (6) 


where 6,,,,, is the Dirac distribution [9] at the point P(z, , z.). An equivalent definition 
is, that for every ¢ which is continuously twice differentiable and zero outside some 
compact subset of A, , 


(Xo , 20) = | / ep] (y) dx dz. (7) 


Thus, if we express ep with the variables £ and ¢ (£ and ¢, are the values of these variables 
for the point P), (7) can be written 


elto , to) = — [fener *M(w)(B/2)r'e™* ae dt 


and é> is a solution of 


M(u) = —(2/3)ro'208¢..¢6 (8) 
similarly, with ~ and X, ep is found to satisfy 
N(u) = (2/8)"*ro'”*8¢,.2. - (9) 


5¢,,:, and 4;,,,, are the Dirac distributions for the two variables &, ¢, and &, \ respectively. 

Although the equation N(u) = 0 does not belong, strictly speaking, to the class 
considered in [7], and [8], it can be studied by the same method. In order to show the 
extension of this method, a simple Dirichlet problem will be considered in the next 
section. 

3. Singular Dirichlet problem for the region x > 0, z > 0. 

A new expression of the Green’s function will now be derived by the same method 
which will be used later for the Tricomi problem. 

In the &, ¢ plane this region is mapped into the strip 0 < ¢ < 1; in the &, \ plane into a 
similar strip 0 < A < A, , (A; being the value of \ which corresponds to t = 0). We 
introduce the Fourier transform U(a, d) of u(é, 4), U = Su, which for summable func- 
tions may be written 


+2 


Sateen / exp (—2imat) ud, u=9"U = / exp (2irat) Uda (10) 


—-@ 


and use, as in [7], the extension of this transform to distributions [9]. It is easy to show 
[7] that the transform of (9) is 
n(U) = Uy + (1 — 0) [(2/8)ixa — 4n°o?JU = (2/8)'*r5'”* exp (2imak,) d;, (11) 


where 6,, is the Dirac distribution of one variable \ at the point \ = A, . In order to 
solve (11), the following notation is introduced: S,(A, a) and S(A, a) are the solutions 
of n(U) = 0 which satisfy 


. te] 0 
Sir , @) = 0, S(O, a) = 0, an SQ, » @) sat 1, an S(O, a) = 1. 
The Fourier transform E> of ep is then defined by Ep = (3/2)7'*ro'”* exp (— 2iwat,) Et 
Et = 0 a)S(o,a), AwM<A<M, 2) 
hS,(ro ; a) S(A, a), 0 < ny < Xo ’ 
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h being such that the jump in the value of the first derivative of E? with respect to d at 


\ = Ap be equal to + 1. Obviously h~’ = S(A,, a) = — S,(0, a). 
On the other hand S, and S are easily found from the first Darboux solutions; set 
8 = 2ima + 1/6, r=, (13) 


one can write 
S,(\, a) = —tF(7/12 + 8/2, 7/12 — 8/2, 3/2, 7), } (14) 
S(A, a) = 3/2(1 — 7)“ F(5/12 + 8/2, 5/12 — s/2, 4/3, 1 — 0), 


where F(a, b, c, r) denotes the usual hypergeometric function. Accordingly h can be 
expressed in terms of gamma functions 


h = 27(11/12 + s/2)T(11/12 — s/2)[T(1/2)T(1/3)}". (15) 


Now we must investigate the most convenient form in which to write the inverse trans- 
form. The equation (11) will be identical to equation (20) of [7] (apart from some obvious 
changes in the notation) if we introduce the variable 8 defined by 


4B? = 4r’a’ — (2/3)ixa (16) 


and choose the branch which reduces for large values of | a | to 8 = a — 1/12”. Thus 
the asymptotic behavior* of the solutions of (11) for large values of | 8 | are given by 
some formulae similar to formulae** (23) and (37) of [7]. The real axis of the a plane 
and the corresponding line Ri{s} = 1/6 of the s plane are mapped into the contour C in 
the 8 plane (Fig. 3). Therefore that the right hand side of (12) is a distribution whose 


‘ 





ee Spe ty 
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inverse transform can be written using an integral in the s plane and therefore the value 


of ep is given by 
1/6+i@ 


ep = (2/3)"(2inrl/?)— [ tra)" ds (17) 


1/6-i@ 


because r = exp £. The reader will recognize in this formula the Mellin transform which 
in fact may be applied to the operator L(u). From such a solution, the expression for 


*The general study of the asymptotic behavior of these solutions is found in [10]. 
**Or the formulae (9) and (10) of [8]. 
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the doublet at the point M,(0, z,) may be derived easily. The Green’s formula for equa- 
tion (1) shows that a solution u defined in x > 0, z > 0 which is zero along Oz, is given by 


uP) = = i oe ne 


0 


This proves that the doublet d,,,(P) in M,(0, z,) is 


re) 
dy,(P) a —21 ax ep(0, 21 ’ To eB (18) 
It is possible to evaluate the right hand side of (17) and (18) by series expansions. 
For instance, if [,(r) is the polynomial of degree p 


$,(7) = F(—p, 11/6 + p, 3/2, 7), 
then 


7-1 


dy,(P) = (2/8)'2to(1 — 10)? D2 40 (3/2 + p)T(11/6 + p)[T(4/3 + p)p!] 


—4/3 —-1\2 2 ° 
r (rita) if 7. Ste, 


xX B,( 70) 


Gr) if So > te. 


In practical applications to boundary value problems the quantity 0/02 dy,(0, 2) is 
very important. The following result is easily obtained 


(3/2)'*41(11/6) (3/2) [e1'(4/3)] rio °F (11/6, 3/2, 4/3, riro’), 


8/8» du.(0, 2) = , if a ee (19) 
(3/2)'/°41(11/6)1'(3/2) [e0'(4/3) ri "73°F (11/6, 3/2, 4/3, rori”), 


if Ts > Te: 


These results can be checked with those obtained by different methods (methods of 
discontinuous integrals or transforms [11], [12], and [13], or the direct method given in 
[2] and [3]. 

4. Green’s function for the Tricomi problem. The same method can be applied 
in order to find an expression for the Tricomi problem of Fig. 2. The domain to be con- 
sidered in the &, ¢ plane is the half plane ¢ > 0 which corresponds to the strip A, > A > 
Az (Az negative) in the &, \ plane. One must again form a solution of (11) with convenient 
boundary conditions. As previously, the solution S, (A, a)—see (14)—is needed. The 
second solution we need in order to define Ep as in (12) is defined by a property of 
asymptotic behavior. According to the general theory, [7, 8] we have two possibilities 
depending on the orientation of the solutions. It was shown that these two solutions, 
which correspond to the two Tricomi problems which can be defined in the domain 
A, , are 1) a function H,(A, a), a solution of n(u) = 0 which is real for 8 = 78’ (6’ positive), 
see (16), and which tends uniformly toward zero for any 4, < A < A, when #’ tends 
towards + ©; 2) a function H,(A, «) which is simply defined by H(A, 8) = H,(A, Be~*"). 
Let us note that (13), and (16) give — 4”°8” = 4”°8” = s’ — 1/36 and s ~ — 229’ 
for @’ sufficiently large. Noting also that H,(A, a) for 6 = 76’ must tend towards zero 
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when A — d, , that is to say t + +, it is evident that H, and H, are given by the 
following expressions, solutions of n(U) = 0 (see [2] p. 11) 

H, = t(r — 177"? *?F(7/12 — 8/2, 11/12 — s/2, 1 — 8, [1 — rJ7'), (20) 

H, = t(r — 1)77"?-*?F(7/12 + 8/2, 11/12 + 8/2,1 +8, [1 — rz]. 
Therefore, if G>’ or Gp’ are the Fourier transforms of gp’ and gp’—these are the two 
functions we are looking for—we may write as in (12), 

Gp? = (8/2) "ro" exp (—2iwak)G¥  (j = 8, 2) 
with 
G* = sa a) H (a » a), A <A<A ’ (21) 
h;S,(\o , OH (Ao , a), Az <= A << Xo ’ 


h; being such that the jump of the first derivative of G¥ with respect to A for X = Ay is 
equal to + 1. Thus h; is the inverse of the Wronskian of S, and H; , then 


h, = 3'?P(11/12 — js/2)T(7/12 — js/2)[2T(1/2)T1. — js) sin w(1/4 — js/2)}"'. (22) 


Now, as in (17), gp’ is given by the Mellin transform 


1/6+1@ 


aS" _ (2 3) (iar) [ G*(rr,")°""”* ds. (23) 
1/6-i@ 

It is not difficult, of course, to verify for this result some properties of the Green’s 
function which have been proved previously [2]. For instance, the ‘“‘symmetry”’ property 
g>’(M) = gy’ (P), and the “inversion” property 

gr '(M) = gp'(r, t3 70 , to) = (ror) gr (ror”*, t3 0 , to). 

For practical purposes, the most important thing is to find the expression for the 
doublet D\;’(P) at the point M, (0, z;) the doublet associated with the Tricomi problems 
in A, . The same argument used previously in order to derive (18) shows that 


~~ ‘Ge 
w.(P) a ax gp (0, 21 » Xo » 20) (24) 
and thus 
- 2 1/3 2 1 1/6+1@ (2) " 
Du (P) = (2) rr’ Qin vn h;H,;(Xo , a) ro ds. (25) 


In many applications, especially in the transonic problems we have in mind, the only 
thing which is important to know is the value of the normal derivative along Oz, a value 
which can be computed with an integral if we know the value of /dz) Dy) (0, 20) = 
K, (2; , 20). This one can be expressed as some kind of generalized hypergeometric function. 
In fact 

P ‘ 
1/4, 1/12, 5/12, im (26) 


2 
K,(2; » Zo) = — (3/2) 4rz7r57”* :f () 
” | 1/4, 7/12, 11/12, 3/4 
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G’,"* being a Meijer’s G-function, (see [14], p. 206). From a straightforward application 


of the formulae given in [14] p. 208, one obtains with the classical notation ,F,(a, ---, 
a, ;b, --- , b, ; z) for the generalized hypergeometric functions. 
First case: 75 < 1; , 

K,(e; , 20) = —(3/2)'*16r0*[3°”"ari]'3F (4/3, 5/3, 1; 5/6, 7/6, rory’). (27) 


Second case: 7, < 1%, 
Kyle; , 20) = —(3/2)'°2[3° ars’? ]"'F2.(7/6, 5/6, 1; 1/3, 2/3; riro”) 


— (2/3)ri?T (1/2) 1 (7/6) [ero’* T(2/3)]"'F (3/2, 7/6; 2/3; riro”) (28) 
— (3/2)'*5ri* [Bars] *.F (11/6, 3/2; 4/3; riro”). 
These formulae seem to be convenient for numerical computations. 
The values of K,(z, , 2.) may be derived from (27) and (28) with 
2oK (2: , Zo.) = 2:Ko(2o , 2). (29) 
Let us conclude this section by recalling that gp’ is the Green’s function for the 
direct Tricomi problem for P in the hyperbolic half plane (Fig. 4), this function has its 


Zz 


} 
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singularities along PP; , PP, (discontinuities proportional to (— z)~'’‘) and its logarith- 


mic singularities along P,P; . Similarly, gp’ is the Green’s function of the conjugate 
Tricomi problem in A, and has its singularities along PP, , PP, , PP; (Fig. 5). The 
doublet Dy,’(P) is the function which allows one to solve the direct Tricomi problem 
where the value of the solution along OC is zero, D\;’(P) the doublet which allows one 
to solve the conjugate Tricomi problem, when the required solution is zero on the 
characteristic at infinity. When r is infinitely small the limiting value of Dy’ (P) (ro 
finite) must be proportional (the factor may depend upon 1,) to the first Darboux solu- 
tion which gives the asymptotic behavior at infinity of the flow around the profile of 
the Mach number 1. This principal value of Dy’ (P) is obtained by considering the 
positive residue in (25). It is proportional to 
To *to(ts — 1)-*°F(4/8, 5/8, 5/2, [1 — t)~’) 

= Kro*{(r — 2)'*82 —r) —¢ +2)'%8z+n}, (30) 


where K is a numerical constant. 
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5. A special generalized Tricomi problem. It was often suggested, for instance 
in [15], that a generalized Tricomi problem could be considered, for which the solution 
is to be found in a domain A, bounded by an arc AB drawn in the elliptic half plane 
with end points A and B on the z axis, an arc AC in the hyperbolic half plane with a 
time-like direction if the direction AB is taken as the time direction, and an arc BC of a 
characteristic. The data for such a problem are the values of the solution u along the 
arc AB and the are AC. By a transformation of the group of the Poincare geometry 
associated with the Tricomi equation, it is always possible to reduce this problem to a 
similar one in which BC is the characteristic at infinity. Now consider the special case 
for which the arcs OB and OC are such that ¢ is constant along each of these arcs. The 
corresponding domain in the £, \ plane is mapped onto a strip A, < A; parallel to the 
t axis. Thus, it is clear that the previous method must allow one to give the solution 
of this generalized Tricomi problem. 

In the &, \ plane we have to start with the fundamental solution e?(M) of N(u) = 0, 
which is “orientated” in the direction opposite to that of the & axis. P(& , Ao) is the 
singular point of e$(M) and the function has a singular behavior along the singular 
characteristics schematically represented on (Fig. 6); e#(/) has a discontinuity on each 
solid characteristic and a logarithmic singularity along each dotted characteristic. On 
the other hand e#(M) will be chosen in such a way that it vanishes when M lies on 
\ = A; orA = A, . The expression of such a solution was given in [7]. Let us call S;(A, a) 
and S,(A, a) the solutions of n(U) = 0 which satisfy the following conditions: 


S3(A3 , a) = 0, Si(\ , a) = 0, 0/dX S3(A3 ’ 0) = E. 0/AX Si, , 0) = :. (31) 
Using this notation, e$(M) is the inverse Fourier transform ( is taken equal to 0 in 
(32)) of 
Et - wees a) S4(Ao ’ a), Xo < A < As ’ (32) 
WS,Q, a) S3(Ao ’ a), N ca < Xo ’ 
with* 
w = S3(A4 , a) = — S,(A3 , a). (33) 


*The functions S;(A, a), S«(A, a) may be easily expressed using the hypergeometric functions of the 
first Darboux solutions. 
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It was shown that expressions (32) are meromophic functions with respect to the variable 
B defined by (16), which have poles for a sequence of real values of 8 and for a sequence 
of purely imaginary values of 8 and that, in order to obtain the fundamental solution 
with the orientation given by Fig. 6, the integral which gives the inverse Fourier trans- 


Ny 








\ Xx / P 
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be ba 








singular characteristics 
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form of (32) must be taken in the 8 plane on a line which leaves the real poles below 
and the imaginary poles with positive ordinates above. Recalling the relations which 
join s and B: s’ = 1/6 — 42°8’, s ~ 2ix® for | 8 | large, it is clear that: 


e#(M) = (2in)™ I Ex(rrs)"-* ds, (34) 
eC, 


where ©, is schematized in the figure 7. It was shown also that e$(M) tends towards 
zero when ¢ tends towards + © for every \, < A <A;. 


é 








Fia. 7, 
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Now, according to (9) the Green’s function gp(M) for the special problem we con- 
sider is 

gr(M) = (2/3)'*ro'“e#(M). (35) 

This function is zero on the characteristic at infinity, on the arc OB and on the arc OC, 

but has a singular point in 0. When P lies in the elliptic part of Ay , gp(M) has the classical 

logarithmic singularity in the vicinity of P, when P lies in the hyperbolic part, gp(M/) 

has a singular behavior along the characteristic lines schematized in Fig. 8. Although 


zh 
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we do not want to enter into all the details, it is quite clear that the previous result 
allows us to state the uniqueness and existence theorems for the solution of this special 
generalized Tricomi problem when the value of the unknown function is prescribed 
along OB(t = t,) and OC(t = t,), by application of the Green’s formula. A small difficulty 
arises in the vicinity of o because of the singularity of the Green’s function near the origin, 
since the application of the Green’s formula gives rise to a series of integrals. But this 
series can be shown to be convergent. Such a situation always arises in classical problems 
for hyperbolic equations in connection with what E. Picard [16] has called a fourth 
boundary value problem, (see for instance [17]). 

6. Concluding remarks. The previous results allow one to simplify a great deal 
the proof of the existence of the solution of the Tricomi problem. To outline very briefly 
such a proof, it is first of all clear that the result of Section 4 gives the direct solution of 
the problem when the contour AMB (Fig. 1) is a normal contour. It is also possible to 
give a direct solution when the contour AMB is, in the z, y plane, any circular arc with 
end-points A and B, by using the result of Section 5. Now the general case when the 
contour AMB is arbitrary can be solved by using the Schwartz process as in [2] and 
the sequence of solutions obtained by this method converges towards the solution by 
application of the maximum principle [2]. 
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DETERMINATION OF COEFFICIENTS OF CAPACITANCE OF 
REGIONS BOUNDED BY COLLINEAR SLITS AND OF RELATED REGIONS* 


BY 
BERNARD EPSTEIN 
University of Pennsylvania 


1. Introduction. In the study of electrostatic field problems the principal objective 
usually is to determine the potential and its gradient (the field strength) throughout a 
given domain bounded by a system of conductors. Frequently, however, it is necessary 
only to determine certain constants of capacitance. In this paper we consider the latter 
problem for a certain class of plane domains. 

Let D be a domain consisting of the entire plane with any finite number of slits 
along a single line. Several formulas for the coefficients of capacitance of such a domain 
are derived, two of which appear to be well suited for numerical computations. One 
of these formulas is based on the explicit representation of the potential as the real 
part of an analytic function [1] while the other formula has the feature of requiring a 
knowledge of the potential only on the line containing the boundary components of D; 
it does not involve any derivatives of the potential. A convenient method for determining 
the field along the line containing the boundary components has been given by the 
author in a previous paper [2]. 

Since the coefficients of capacitance of a domain are invariant under conformal 
mapping, the formulas which are derived may be employed to compute the coefficients 
of any domain which can be conformally mapped upon a domain D of the type described 
above. This procedure can be applied to a certain class of domains which are of practical 
interest. In these cases the mapping problem involves essentially only simply-connected 
domains rather than multiply-connected domains. One particular case of interest, the 
‘bi-filar shielded cable’ is considered in some detail, and as an illustration of the procedure, 
the coefficients of capacitance are evaluated numerically for one such domain. 

2. Regions bounded by collinear slits. We consider here the problem of determining 
effectively the so-called coefficients of capacitance of a domain D consisting of the 
extended (z, y)-plane with a finite number, m, of collinear slits, cut along what may be 
assumed to be the z-axis. 

We recall the definition of these coefficients: 


Ou; , 
as de, (j,k = 1, 2, ~~ , mM), (2.1) 


Ce 
where u,; , the harmonic measure of the jth boundary component, is the harmonic function 
whose boundary values are unity on the jth component and zero on all other components 
(the existence and uniqueness of this harmonic function is well known from the theory 
of the Dirichlet problem); C, is any curve, described in the positive sense, surrounding 
only the kth component; and @/dn indicates differentiation in the direction of the outward 
normal. As follows immediately from the Cauchy-Riemann equations, the p;, may be 
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defined alternatively as the increment in the harmonic function v; conjugate to u; which 
results when C, is described once-in the negative sense. 
We also recall several important properties of these coefficients of capacitance: 


Die = Dai ; (2.2a) 

> Pn = 0, j=1,2, +++ ,m; (2.2b) 
Di; > 0; (2.2c) 

Pix <0, JGH#k. (2.2d) 


The p;, are, of course, conformal invariants. Hence we may assume that D consists of 
the entire plane minus a finite number of slits lying on the z-axis, one of which extends 
to infinity in both directions (see Fig. 1); this configuration can always be realized by 











Fie. 1 


a suitable inversion. It will be seen that this assumption eliminates the possibility of 
any convergence difficulties. For brevity such domains will be called slit-domains. We 
number the finite slits 1, 2, --- , m — 1 from left to right; and the infinite slit is the mth. 

We proceed to derive various formulas for the quantities p;, which may prove useful 
for numerical computations. First we shall employ the definition (2.1); later we shall 
derive a formula based on the alternative definition given following (2.1). 

Taking into account (2.2a) and (2.2b) we see that it suffices to determine the p;, 
with 7 < m, k < m. Let f;({x) = u;(z, 0), so that, in particular, f;(z) = Oforz < a,, 
x > b,, . By the Poisson formula we have, for any point off the z-axis: 
ly| wf ) __ fi) dé a (2.3) 

. — x)” +y 
We take as the curve C, a rectangle with vertical sides passing through the gaps (a, , b,) 
and (a,.; , b4+:). It is easily found from (2.3), by differentiating under the integral sign, 


that for any value of z satisfying a, < z < D,, 
| | 9 
Ou; | 2A“ | du; | 2A 
—|<— —|<-5 A= x (| a, |, | 0, | 2.4 
az |—alyP’ lay | Say?’ { max (| a, |, | bn |) (2.4) 


u;(z, y) = 


(one uses, of course, the fact that | f;(€) | < 1). Hence as the horizontal sides recede to 
infinity the contribution to the integral (2.1) from these sides vanishes, and therefore 
(2.1) may be written as follows: 
Pik = Tin ick | ve ’ (2.5) 
where 
* Ou; 


I,,= on WY a,<2 <b, (2.6) 


(of course pin = I;,m — I;,1). 
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We wish to rewrite J;,, in such a form as to involve only the values of u; on the 
x-axis, i.e., the values f;(x). If one expresses (du;)/dx by differentiating (2.3), inserts 
this expression into (2.6), and interchanges the order of integration, one obtains the 


following expression for J;,, : 


2p [-LOe * <28<i&. (2.7) 


Here the symbol P denotes the Cauchy principal value of the integral. 
A second expression for J;,, is cbtained by the following artifice. Since the right side 
of (2.6) gives the same value for any choice of z in the rth gap, we integrate both sides 
f (2.6) over this gap. Then we obtain 


(b, — a,)1;, = Fr St 7 dx. 


Interchanging the order of integration [this is ae justified with the aid of the first 
inequality in (2.4)] we obtain 


* 2 aw 
(,-a).= [| Sut da dy = | {u(b, ,y) — ula, , »)} dy 


Tw ™ 1 _ 1 
- [. ly | {f £0| z OR, b,)? + y (é oe, a,)” + = | as} dy. 


Now we interchange the order of integration once more, and obtain the following ex- 


(2.8) 








Il 


pression for J,,, 





lin ao I f@) mn | E= (2.9) 


A third expression for J;,, is obtained as follows. Since -" sai is continuous and f;(é) 
and In | ( — a,)/(€ — b,) | are absolutely integrable, an integration by parts enables us 
to rewrite (2.9) in the form: 


lie = Ze ay | HOWE - a) In | - 


a(b 


—(€—b,) n|g—}, |] de (2.10) 





the integrated term vanishes since f;(a,) = f;(b,) = 

Equations (2.7), (2.9), (2.10), together with (2.5), give three formulas for the co- 
efficients p;, which employ only the values of u; on the z-axis. These equations were all 
derived using the definition (2.1) of p;, . A fourth formula is obtained by using the alter- 
native definition, as follows. It is shown in [1], Sec. 91 that the analytic function 
w;(z) = u; + tw; must satisfy the condition 


P (2) (2.11) 


wiz) = -—— =) 
‘- Il@-ae- i} 


r=l 





where P;(z) is a polynomial, of degree not exceeding m — 2, with m — 1 real coefficients 
which are uniquely determined by the conditions*: 


+1, for r=j, 
be 
| wi(z)dz=4-1, for r=j+1, (2.12) 
0, for r#¥j, r#jt+il. 


*It must be remembered that the denominator of w;(z) changes sign on successive intervals (a, , br). 
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(It is shown in reference [1] that of the m conditions (2.12) only m — 1 are independent.) 
Now employing the alternative definition of p;, and taking for C, the doubly-counted 
kth interval, one easily obtains from (2.11) the result 


Da = ¥2 [~ —< (2.13) 
. { IT & — a)¢ - i} 





where the ambiguity in sign is most easily resolved with the aid of (2.2c) and (2.2d). 

Insofar as numerical computations are concerned, it would appear that formulas 
(2.9) and (2.13) are especially suitable. While (2.13) has an advantage over (2.9) in 
that integration is necessary only over a single interval, it involves the solution of the 
system (2.12) for the coefficients of the polynomial P;(z) and this may become laborious 
for large values of m. In this case it might prove preferable to employ (2.9), obtaining 
the function f;(z) to the desired degree of approximation by the method given in [2], 
Sec. 6. 

3. Extension to more general regions. Since the coefficients of capacitance are 
conformal invariants, they may be determined by the method described above for any 
domain which can be mapped conformally onto a slit-domain. A simple example of 
such a domain is the entire (z, y)-plane slit along a finite number of arcs of a circle (see 
[1], Sec. 92). By a suitable linear transformation we can map the circumference of this 
circle onto the z-axis, thus obtaining a slit-domain. 

Another domain to which this method applies is a domain bounded externally by 
one circle and internally by two others. By a suitable linear transformation such a 
region may always be mapped into a domain with the centers of all three circles collinear 


(see Fig. 2). 








Fic. 2. 


Now suppose that the upper half of this domain, which is simply connected, is mapped 
conformally by an analytic function ¢ = f(z) upon the upper half of the ¢-plane. Then 
the upper halves of the three circles are mapped into segments of the real axis of the 
f-plane, and the three segments of the axis of symmetry (AB, CD, EF) are mapped 
into the remainder of the real axis. By the Schwarz reflection principle, the entire domain 
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of Fig. 2 is then conformally mapped onto the {-plane with three slits along the real 
axis, i.e., upon a slit-domain. A particular case of practical interest in the ‘bi-filar shielded 
conductor’. This case will be discussed in some detail in Sec. 4. 

Closely related to the configuration of Fig. 2 is that of the plane with any finite 
number of circular apertures the centers of which all lie on one line (see Fig. 3). Each 


y 





Pe 
VY 


O0-o 


Fia. 3. 


ax 
a, 





half of this region is simply connected. If, as before, we can make a conformal mapping 
of the upper half of this region upon the upper half of the ¢-plane, then, as in the previous 
case, the reflection principle immediately gives the mapping of the entire configuration 
upon the w-plane with collinear slits. 

The use of the Schwarz reflection principle to transform a region into one bounded 
by collinear slits may prove of use in studying the electrostatic field created by an 
electrified grid. 

4. An illustrative example. In the case of the domains described in Sec. 3 it has 
been shown that the problem of determining the coefficients p;, may be reduced essentially 
to the problem of finding a mapping function for a simply connected domain, regardless 
of the connectivity of the original domain. Since the simply connected domain in question 
is bounded by parts of circles and straight lines, the required mapping can be obtained, 
in principle, by solving a certain non-linear third-order differential equation which is 
closely related to a certain linear second order equation (see [3], pp. 198-208). However, 
in practice this fact is of little use, for these differential equations contain the so-called 
accessory parameters of the domain, which cannot readily be determined except in the 
most elementary cases. But for regions of the type under consideration, it might be 
possible to overcome this difficulty with the help of some of the many approximate 
conformal mapping techniques that have been developed. 

Here, by way of an example, we shall consider a case of the ‘bi-filar shielded cable’ 
which can be treated by elementary methods. 

Referring to Fig. 4, let the radius of the outer circle be taken equal to unity (this 
assumption does not reduce the generality), and let the points B, C, D, E correspond 
toz = —£8, —a, a, B respectively. By the mapping function 


t= (as)"(z + 1) (4.1) 


the interior of the outer circle is mapped upon the exterior of the slit —2(ag)'’” < 
t < 2(a@)'”’, and hence the domain bounded by the three circles is mapped onto the 
¢-plane minus the aforementioned slit and the images of the two interior circles. In 
particular, if 8 < 1, so that the inner circles are small and located close to the center of 
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the outer circle, then the mapping (4.1) may be approximated in the neighborhood of 
the inner circles by the linear mapping 


1/2 
p= (4.2) 
oe 
P f nt, 
/ | i 
OO 
\ | / 
\ | F / 
Fic. 4. 


which transforms the inner circles into circles which meet the real axis at +(a/g)'” 
and +(8/a)'”. Hence the domain of Fig. 4 is mapped by (4.1) approximately into 
the region indicated in Fig. 5. The centers of the circles in Fig. 5 are given by ¢ = ta, 


C-plane 


“par Jay -2(apy"* 208)” ae 


Fia. 5. 








where a = (a + £)/[2(a8)'”], and the radius of each circle is given by R = (8 — a)/ 
[2(af)'’”]. Since a? — R? = 1, the exterior of these two circles may be mapped upon the 
exterior of two slits of the real axis by the function (see [3], p. 297, example 9) 


_ 1+ fli + 9/0 — 9, oe] 
~T=fla-+ 9/0 — 0), ol aa 





Ww 


where 
a (8)? ae (a)*/ 
(6)* + @) 
and 


‘ Il (1 4+ pr” *2’) Il (1 + a 
f(t, ) =" = a2! (4.4) 


t o o 
I] (1 + gp) I] (1 f- p**e-) 
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(see [3], Chap. VI, Sec. 3, Eqs. (32) and (49)). The circles map onto the twice-counted 
slits 


_1+ Lp) _1— Lp), 
1— Lip) =” = TF Lp)’ 


and 

1 — L(p) 1 + L(p) 

coe Con fo Skee 

1+ L() =" =1-Le) 
respectively, where 

@ 1 + p” ) 
= 
L(p) = 2p I] [toe ; (4.5) 


The interval —2(ag)'”? < ¢ < 2(a@)'” maps onto the infinite interval 


1 — f[{1 + 2(ap)'”}/{1 — 2(aB)"”*}, p] 
1+ f[{1 + 2(aB)'*}/}1 — 2(ap)'”*}, p] 


of the real axis. Thus the domain of Fig. 4 has been approximately mapped into a slit- 
domain, and the formulas given there may be employed to determine the p;, . 

To give a numerical! illustration, let a = 1/20, 8 = 1/5. Then p = 1/3, and, since 
the infinite products (4.4) and (4.5) converge rapidly, one obtains very easily 





|w|> 








1 — L(p) 1 + L(p) 

———e ws ().212, —————~ ws 4.726, 

1 + L(p) 1 — L(p) (4.6) 
es fifi + 2(a)""*}/{1 cox 2(a8)'”?} , p} = 8.628 
1 + f[{l + 2(@s)'*}/{1 — 2(@s)'*}, p]) 

In order to work with more convenient numbers, we employ the transformation 
2 
a oe (4.7) 


thus obtaining in the Z-plane a slit-domain characterized by the following numbers 
(see Fig. 1): 


a; = —8.1 b, = — 4.5, a» = —0.2, b, = 0.2, a3 = 4.5, bs = 8.1. (4.8) 


] 


On account of the low connectivity of the domain (m = 3) it was decided to employ 
formula (2.13) to evaluate the p;, . The coefficients p.. and p,, (the numbering of the 
boundary components is given in Fig. 4 and Fig. 5, and is in accordance with the number- 
ing given preceding Fig. 1) were determined as follows. The polynomial P,(Z), which 
according to the discussion of Sec. 2 is of first degree, was replaced in (2.12) by the 
quadratic polynomial aZ* + bZ + c with unknown coefficients a, b, c. The integrals 
appearing in (2.12) were computed by interpolation (the Tchebyshev five-point formula 
was used) and the resulting system of linear equations was solved for the coefficients 
a, b, c, with the following results: 


a=-93X10°, b= —1.962, c= —11.608. (4.9) 
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(The radical appearing in (2.12) was taken positive on the first and third intervals and 
negative on the second interval.) The extremely small value obtained for a, whose exact 
value must be zero, serves as an excelient check on the accuracy of the computations. 
The expression for P,(Z) thus obtained was then employed in (2.13), with the following 


results: 
Pan = —2.153, Poo = +3.820. (4.10) 


By symmetry, p;; = P22 , and now the remaining six coefficients are obtained with the 
aid of (2.2a) and (2.2b). Thus, the following table of values of the p;, was obtained. 














\ J 

k \ 1 2 3 
1 |  +8.820 | ~2.153 | 1.667 
2 ~2.153 | +3. 820 —1.667 
3 | —1.667 | —1.667 +3.334 





As a check, p.. was also computed by means of (2.5) and (2.10), and excellent agreement 


was obtained. 
Acknowledgement. We wish to thank Mr. Charles Kahane, who so capably carried 


out the computations connected with the illustrative example of Sec. 4. 
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GENERAL THEORY OF ELASTIC STABILITY* 


BY 
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Harvard University 


Summary. Some general topics in elastic stability are discussed. In particular, 
attention is given to the relationship between adjacent-equilibrium-position and energy 
techniques, to the effects of non-linearity, and to the sensitivity of certain stability 
problems to the character of the loading. 

1. Introduction. In analyzing the stability of an equilibrium state of a particular 
elastic system, those terms which arise from the equilibrium condition eventually cancel; 
consequently a number of writers have found it desirable to discuss stability in a general 
manner, removing these terms once and for all, and directing their attention towards 
the remaining terms. Usually, a certain equilibrium state is postulated, and one of two 
criteria is then used to determine the stability of this state. The first criterion states that a 
structure is unstable if an adjacent equilibrium state exists, whereas the second requires 
for instability that the over-all potential energy not be a relative minimum. 

In setting up these criteria in analytical form, it is recognized that some sort of 
non-linearity is essential, in order for example to evade the uniqueness theorem of linear 
elasticity. Such non-linearity may arise either from the geometry of the situation (large 
displacements or non-linearized boundary conditions) or from the inclusion of higher 
order terms in the stress-strain law. Although the above uniqueness theorem reason 
may not be valid (on the grounds that the usual proof of this theorem contemplates the 
same body configuration for the two supposedly different stress-strain states to be 
proven identical, and so is not applicable to stability problems in any event), it is never- 
theless clear that because of the cancellation of equilibrium terms it is well to include all 
important higher order terms. The method of incorporating these terms varies widely, 
as will be seen from a study of treatments by Bryan [1], Southwell [2], Biezeno and 
Hencky [3], Trefftz [4], Biot [5], Neuber [6], Prager [7], Goodier and Plass [8], and others. 
Most of the assumptions concerning non-linearity made in these treatments seem rather 
artificial—for example, Trefftz and Goodier obtain non-linearities by regarding as funda- 
mental a curvilinear coordinate system which moves with the material fibers of the body, 
and there seems to be little basis for this method. Similarly, Prager uses techniques 
of superposition which are questionable when dealing with non-linear effects.** 

The role of non-linearity has been further complicated as a result of a paper by 
Goodier [9], in which it is stated that the correct equations for the torsional buckling of 
a bar can apparently not be obtained by conventional energy techniques. Goodier gives 
a rather complicated analysis of this problem, using a Trefftz-type method, and his 
incorporation of non-linearities again appears to be quite arbitrary. In the problem of 
shell buckling, as discussed for example by Karmdn and Tsien [10], it has been suggested 





*Received May 25, 1955; revised manuscript received September 22, 1955. 

**Nevertheless, the treatment of Prager is probably the clearest available. He uses the adjacent 
equilibrium position technique of Biezeno and Hencky, but obtains their results in a much more compact 
manner. In addition, the effects of inelasticity and thermal gradients are considered, and the final eigen- 
value problem is thrown into the form of a variational principle. 
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that some sort of non-linearity may be partly responsible for the discrepancy between 
experiment and theory (besides the known effect of initial irregularity). 

On the other hand, it is known that the effect of non-linearity in the stress-strain 
laws governing the usual structural metals is only of the same order of magnitude as the 
uncertainties in the ordinary elastic constants, and it does not seem physically reasonable 
that such effects should materially influence practical stability problems. In the present 
analysis, the stability problem is first analyzed without approximation; an engineering 
approximation is then obtained by consistently neglecting terms of a certain order of 
magnitude. Specifically, an arbitrary elastic body in equilibrium under certain loading 
is considered. An arbitrary virtual displacement is assumed, and the discrepancy between 
the work done by the loading and the increase in internal energy is calculated by means 
of the exact stress-strain relationship of non-linear elasticity; the most convenient 
form of this law is that due to Murnaghan [11]. It is decided that a positive discrepancy 
is the necessary and sufficient condition for instability, and the analytical consequences of 
this are worked out. The appropriate engineering approximation is made in the final 
result, and it is found that the result is different from and simpler than those usually 
obtained. As a special case, the problem of Goodier [9] is considered and it is found 
that the correct equations are obtained in a straightforward manner. 

Recent discussion by Pfliiger [12] and Ziegler [13] have directed attention towards 
certain fundamental problems in stability. Ziegler has shown by exemplification with 
non-conservative systems that the result of examining for stability the equations of 
motion of a perturbed mechanical system do not necessarily coincide with the familiar 
energy or adjacent-equilibrium-position techniques. Since the equations-of-motion 
method must be regarded as basie, Ziegler’s work gives rise to some doubts as to the 
usefulness of the other methods. However, for the case of conservative systems (to 
which we restrict ourselves for the present; plastic buckling will be discussed elsewhere), 
the equations-of-motion method and the energy method are equivalent (a proof will be 
found in Whittaker [14]), and so the energy method may be used with confidence. How- 
ever, the adjacent-equilibrium method and the energy method are certainly not equivalent 
even for conservative systems. Consider for example an (always elastic) column com- 
pressed beyond the buckling load but restrained from buckling; if the constraints are 
removed the column will buckle despite the absence of an adjacent equilibrium position. 

But we can perhaps obtain an equivalence by altering the problem somewhat. 
Consider an elastic system which is in stable equilibrium under certain loading. As the 
loading is increased in some manner, the system following an equilibrium path, a point 
of instability (by the energy criterion) may be reached; does an adjacent equilibrium 
position exist at this critical point? This is a question of considerable practical signifi- 
cance. For example, in the problem of “tin-canning’’—1i.e., the problem of stability of 
a fairly flat shell under lateral pressure, where at a certain critical load the shell tends 
to snap suddenly through into an entirely different equilibrium position—there seems 
to be no adjacent equilibrium position (in the conventional Euler-column sense) cor- 
responding to the critical load. Can an adjacent-equilibrium-technique then (as used in 
practice) give the correct answer to such problems? 

It will be shown that the two techniques are indeed identical for the altered problem 
of the last paragraph (so that, for example, the correct answer in tin-canning problems 
is that conventionally obtained, and one must look elsewhere for the discrepancy between 
theory and experiment). In particular, it will be found that the differential equations of 
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first variation of the general energy principle are precisely the same as the conditions 
for existence of an adjacent equilibrium position, these conditions being again calcu- 
lated without approximation. 

Another stability topic of considerable interest relates to the precise character of 
the loading applied to an elastic system. Such effects are not considered in previous 
general stability analyses, yet have been shown to be important by Tsien [15]. A particu- 
larly subtle example is given here in which, for any perturbation, the first-order work done 
by two alternative types of loading is the same, yet the buckling loads are widely different. 
Since the problem is a very practical one (buckling of a long cylinder under external 
pressure), it is clear that the character of the loading must be included in any general 
stability criterion. The appropriate analysis will be given for the two types of loading 
of greatest importance, viz., dead loading and pressure loading. 

2. Analytical condition or instability with dead loading. Consider an arbitrary 
elastic body which is initially free from stress (state I). By the application of load or of 
heat, the body alters position and shape (and achieves state II). The material particle 
initially at the point (a, , a2 , a3) has now moved to (2, , Zz , 23), where subscripts are 
used to distinguish between the usual three fixed Cartesian axes. The ith component 


of the displacement vector is given by 


9 = 2% — Ga; (1) 


and the Lagrangian strain tensor is 


Nii = 11 dv, ‘Oa; + dv; /da; + (dv,/da;)(dv,/da;)], (2) 


where the summation convention is used (here and in the future) for repeated subscripts. 

If U is the internal energy per unit mass, a symmetric function of the nine ,; , of 
the absolute temperature 7 (or entropy S), and of position, then by Murnaghan’s 
treatment [11] the Eulerian stress 7,;; in state II is given by 


i ies p(d U/8np) (02; /da,) (dx;/da,), (3) 


where p is the density in state II and the partial differentiation of U is to be carried 
out at constant entropy. It is now required to analyze the stability of the body in its 
deformed state II. From Sec. 1, the body will be considered stable if for each infinitesimal 
displacement (compatible with the Loundary conditions) the work that would be done 
by the surface and body forces does not exceed that absorbed as an increase in internal 
energy.* If this condition is not met, then for some virtual displacement excess energy 
would be available for use as kinetic energy, and the appropriate displacement will 
increase in magnitude. 

The body force per unit mass, F; , will be assumed constant (e.g., gravitation). The 
surface loading, 7’; per unit area in state II, is considered to be produced by fixed loads 
which vary neither in total magnitude nor direction during the trial displacement. 
Thus, under such “dead” loading, the material particles constituting a portion of the 
surface in state II will always be subject to the same total surface vector force, irrespective 
of their orientation or total area, throughout the trial displacement. Consequently, the 





*In the equivalent potential energy form, this is the usual energetic stability criterion. We use the 
above form because of its additional generality; as will be shown elsewhere, it can then be applied to 


certain non-conservative systems also. 








136 CARL E. PEARSON [Vol. XIV, No. 2 


work done by the body and surface forces in a displacement u; from state II would be, 
exactly, 


w < [ Fu. av + / Tu; dS, 


where the volume and surface integrals are calculated for state II. Altering to volume 
integrals and using the equations of equilibrium gives 


vs / ,,(du,/dz,) dV 
which, upon substitution from Eq. (3), becomes 
W = | @u/an,..(az. /aa,)(du,/da,)p aV. (4) 
The increase in internal energy is, exactly, 
a = [ (U’ — U)pav, (5) 


where U’ denotes the internal energy per unit mass following the displacement u,; , 
and depends on the temperature of that state as well as on u,; . Note that the volume 
integral is still calculated for state II (this is allowable because the element of mass, 
p aV, is invariant). 

Consequently, the general condition for stability is that, for each allowable u, , 


e 


| ((8U/8n,.).(02,/Aa,)(du,/da,) — {U’ — U}]pdV <0 (6) 
(in particular the integral vanishes for u; = 0; if it vanishes for some u; ~ 0 but is 
non-positive for all u,; , the state will be called neutrally stable.) 

It is now necessary to calculate U’. Because buckling is usually rapid, it is reasonable 
to require the displacement u; to be of an adiabatic character, and we will make this 
assumption. Had we at this point insisted on an isothermal motion, an entirely analogous 
calculation (best carried out by use of the Helmholtz free energy function instead of U) 
could have been made, and the same final results would be obtained in the sequel except 
that the isothermal rather than the adiabatic elastic constants would appear. Experi- 
mentally, the difference between these constants is negligible; then, using the fact that 
in general the motion u; of the body would be somewhere between adiabatic and isother- 
mal, it is seen that the particular thermal assumption at this point makes little difference. 
In any event, we consider for definiteness an adiabatic motion, so that in the power 
series expansion of U, viz., 

U’ — U = (8U/Am¢) SMpq + 2(8°U/Ay¢ 90:5) 50:5 SMe + *°- (7) 
all partial derivatives are to be calculated for constant entropy and for state II. Using 
5n:; = 3[(0x,/da;)(du,/dat) + (dx,/da;)(du,/da;) + (du,/da;)(du,/da;) } (8) 


in Eq. (7), and substituting the result into Eq. (6), gives as the condition for stability that 


/ p dV(8U/an,;)(du,/8a,)(du,/8a,) + (8°U/Ang; On.) Onc; Btpe + °°°] 
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be greater than zero for all non-zero permissible u; . Alternatively, use of Eq. (3) allows 
the condition to be written as 


[ aV {r5e(du,/8x,)(du,/2x,) + o(8°U/dn;; Opa) 5ni; ONp¢ + °° *] > 0. (9) 


Except for pathological cases, only second-order terms need be considered, and the 
criterion for stability becomes 


/ dV[r,,(du,/dzx,)(du,/dzx,) (10) 


+ p(d°U/dn;; On,,)(dx,/da;)(Ax,/da,)(du,/da;)(du,/da,)] > 0. 
Here, all quantities are calculated for state II. For adiabatic virtual displacements, 
this criterion is exact, and must be used wherever non-linearity of the stress-strain law 
is essential. 

3. Engineering approximation. Isotropic media. The second term in Eq. (10) 
may be calculated by means of a power series expansion in 7;; in terms of the various 
derivatives of U evaluated for state I. Using the subscript “0”’ to indicate state I, 

a°U Oni; ONve - (0°U/dn,; ONe)o + (8°U/dn;; ONpa On, s)oNrs + etea (11) 


For structural metals, the magnitude of the second term in Eq. (11) is generally smaller 
than the uncertainty in the experimental value of the first term (the first term repre- 
sents the usual elastic constants, and the second and following terms represent non-linear 
elastic effects). Consequently, it is reasonable to approximate the second term coefficient 


by 
Po (8°U/dn;; ONpa)o ’ 


where the density in state II has also been replaced by the density in state I. This term 
is recognized as the conventional (adiabatic) elastic coefficient and will be denoted by 
Ciipg » Then the second term may be written 

[c? n¢(0x,/0a;)(x,/da;)(Ax,/da,)(O2_,/Aa,) |(du,/Ax,)(Ou,/dX,). (12) 
Now the deformation (although not the displacement) between states I and II is assumed 
small; this means that the partial derivatives inside the square bracket of (12) represent, 
within the approximation being made, a pure rotation. But the quantities c?;,. form a 
Cartesian tensor, so that the quantity in square brackets reduces simply to the elastic 
coefficients for the orientation of state II, i.e., to c,.._, . Thus the stability condition, 
Eq. (10), becomes 


/ dV [r,<(8u,/82,)(8u,/82,) + Creamlrom] > 0, (13) 


where e;; = 3[(du;/dz,;) + (du,;/dx,;)] and where the symmetry property of c,,,,, has been 
used. For isotropic media, 


9 
Crtem _ ol 8.8.0 + 5105 em + —s in | (14) 
1 — 2¢ 


where G is the shear modulus and oa is Poisson’s ratio. Use of this relation gives the 
stability criterion as 


/ av| r(Ou,/22,)(Ou,/82,) + 26 ener + —" cutan} | > 0. (15) 
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The only terms which have been neglected in the derivation of Eq. (15) are those 
which are of higher order in powers of 7;; . 

4. Euler column. Before proceeding with the general theory, it is worthwhile to 
consider a simple example. Let a slender column of length L and cross-sectional area A 
be placed so that its neutral axis coincides with the x; axis, and so that it may buckle 
in the x, — 2; plane only, the ends being restrained from lateral motion. A total load 
P is applied to the end of the column, producing a stress of 73, = —(P/A), all other 
7;; = 0. We now use Eq. (15), and see how large P must be, for certain trial displace- 
ments, before the left-hand side of Eq. (15) becomes negative; such a situation would 
correspond to buckling. Since the trial displacement will usually not be the exactly best 
ones for this purpose, the buckling load obtained in this manner will always be too high; 
this remark clearly holds in general also and is not restricted to the Euler column case 
(see Ref. [8]). In fact, the buckling load P, is given by 


: ~  f fe.m@em + (0/1 — Zoe: :emm} AV . 
P./2GA) = min ——— : - — . 16 
(P./2G . J (Ou, /0x3)(Ou,/dx3) dl (16) 


{ug} 





Whenever P exceeds P, , Eq. (15) shows that the column is unstable. 
Let us choose a trial displacement, being guided in our choice by the tendency of 
’ 5 A 
plane cross sections to remain plane and perpendicular to the neutral axis during bending: 


Uy = U3), 
in = O, (17) 
Us = —2,U'(25), 


where u(2;) is an arbitrary function of x; which vanishes at x, = 0, L, and where u’(z3) 
is its derivative. A straightforward calculation using Eq. (15) gives that 


aL oL as L 
—(P a) 4 | (uj? +1 | ws | + 2¢( == \r i (u’’)? <0 (18) 


for instability, where J is the appropriate moment of inertia of the cross section. Because 
(P/A) < G, the second term in the square bracket is omitted and we obtain 


> 


* wvl-—e ; hy (u’’) 

P. = 2GI\————] min — ’ 
l — 20 (e (u’)? 

(19) 


l—a/) 
= 2G1| —— ( 4) 
= 1(; — ao) w/t 


This answer is too high, because 


so that the displacement (17) is deficient in some respect. The deficiency lies in the fact 
that u, and u, do not contain terms allowing for lateral expansion of the column during 
bending. Actually, instead of setting e,, = €2. = 0 in Eq. (15), we should more correctly 
have set €;, = €22: = —oé3; . A simple calculation shows that the incorporation of such 
terms does not materially alter the first term in Eq. (15). A suitable altered displace- 
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ment would in fact be 
o 2 2 
“wM=uts (xz; — 22)u’’, 
(20) 
Use = o2,2u"', 
le = —2,uU’. 


If this displacement is inserted into Eq. (15), and the magnitude of various terms ex- 
amined (which is most easily done by assuming that u is not far removed from that u 
used to minimize Eq. (19), viz., sin (rz/L)), it is found that the first term of Eq. (15) 
is essentially unaltered, whereas the second becomes (very closely) Ee;; . Then a similar 
salculation to that of Eq. (19) gives 


P, = EI(x/L)’. 


It will be noted in Eq. (15) that the second term is the familiar strain energy term, so 
that the first term must in a sense represent work done by the loading. It is therefore 
not surprising that minor transverse alterations in the u; (these alterations incidentally 
vanishing on the neutral axis) do not affect to any extent the value of the first term. 
This is a rather useful point to note, because it means that simple displacements of the 
type (17) may be used in many column problems, provided only that e,,; and eé22 are 
set equal to (—cé33) in Eq. (15). 

5. Flexural buckling. Goodier [9] has examined the use of energy techniques in 
the flexural buckling of a twisted bar. Since his analysis is geometrically complicated 
and physically questionable it is worthwhile to show that the correct equations are 
obtained by use of Eq. (15) in a routine manner. Since the only question is as to whether 
or not certain terms occur, it is only necessary to consider a simple special case—that 
of a circular cylinder. If the central line coincides with the x,;-axis and if the angle of 
twist per unit length is 0, the stresses are 


3. = —G6x, ’ Ts) = Gx, ° 

Assume a displacement of the form 

u,; = u— Br, 

U2 =v+ pr, 

U3 == —2z,u’ _ 0, 
where u, v, 8 are functions of x; . (If the cylinder were non-circular, a term involving 
the warping function multiplied by 6’ should be adjoined to u, .) Then using the tech- 
nique of Sec. 4, condition (15) requires for stability 


L L 


2G6 | (—w’v''I, + v’'u’'I.) + GIy 


L 
(e’)? + E / Taw’)? + T10")"] > 0 


and minimizing* this expression yields 


EIu’"’ + Mo’ =0,  Elw’’ — Mw’ =0, 8” = 0. 


*That minimization is the appropriate procedure will be shown subsequently. 
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With the appropriate boundary conditions, these coincide with the final results of Goodier. 
(Note that 7, and J, are defined in an opposite way to that of Goodier. Here J, is defined 


as being about the z,-axis, i.e., f 22 dA.) 

6. Curvilinear coordinates. Buckling of a cylinder under dead load. Very often, 
the appropriate coordinate system is not Cartesian; in such cases it is useful to have 
available a more general formulation of Eq. (15). Let the differential element of distance 
be given by 

ds” = h? dy? + h? dy? + h? dy3 , 
where h, , h. , h; are functions of the three curvilinear coordinates y, , y2 , y3 . Then 
denoting by 7,; the curvilinear stress components and by wu; the curvilinear virtual 
displacement component (i.e., in the parametric direction of y;), the first term of Eq. 
(15) may be shown by direct calculation to become 


T {u, uu 
bm | h pura + <r h,, wie 
_ te 1 h; (21) 


UU 2u 2u Qu Ul 
— (i= * & » t+ — i —_ - a a 7 } 2 hy cha») |, 
( ly Umlla 


where a comma indicates iitciabetien with respect to the appropriate y,—thus u,., 


means (0z,/dy,). 
The form of the second term is unaltered, but e;; must now be interpreted as a cur- 


vilinear strain component, perhaps most conveniently given by 


1 h; o h; o (us) i Us $s | - 
-}5[h3 ay, (x Th. oy, \ad * 79 th, oy, t - 


Consider for example a long thin circular shell of mean radius R, under the action 
of an external pressure P of the present dead-loading (see Sec. 3) type. As cylindrical 
coordinates, choose y, = x along the axis of the shell, y, = 0, the polar angle, and y; = r, 
the radial distance from the central axis. Then 


ds’ = dy; + y3 dy2 + dy; . 





Choose 
u, = 0, 
Ul, = 0 += (v — w’), 
R 
uz; = UW, 


where v, w are functions of @ representing the motion of the central surface of the shell, 
and z = r — R. From shell theory, it is known that displacements of this type are suit- 
able for calculating all strains except é€;3 , €23 , 33 - The former two are conventionally 
negligible, and the latter is usually calculated by assuming the induced 73; stresses to 
be much smaller than the bending stresses in the shell. Then the e;; to be used in Eq. 





(15) are 
Cir = Cig = C13 = C23 = O, 
en = E+ 2 pw tw"), 
oC 
€33 = a €22 ’ 


1956] GENERAL THEORY OF ELASTIC STABILITY 141 


where a minor approximation has been made. Using Eq. (21), the stability condition 
(15) becomes that 


| av{(—PE) Cu. + us)? + (te — ta,s)"] + 


tr° 





E 
= o a} > 0, 


where ¢ is the thickness of the shell. Because (PR)/t « E, the first term may be submerged 
in the last to give, approximately, 


, y ’ z — anl\]2 E v’ wow fF ” F 
[aie - wy +50 wr +25[f+E Zw +0") |} > 0. 





Taking a unit length of cylinder and integrating over the cross-sectional area gives 





a E —" t {o" (v’ + w)? + (#/12R’) [5" (w+ w’’)? 
(l-—o)R 2" (v — w’)’ 
which is obtained (very closely) by setting v’ = —w and v = sin 20(v = sin @ would 


correspond to rigid body motion). Then 


————— 
Fe oa (1 oN o)R° ’ (24) 
where J = ¢*/12. Since the usual result is (3/4) of this, it is clear that the assumption of 
dead loading has materially altered the critical load. Load-type sensitivity has been 
remarked for this problem by Stevens [16] and more generally by Tsien [15]; we consider 
it here to exemplify the manner in which the stability criterion will be generalized. 

7. Pressure loading. We return now to the general theory of Sec. 2, and examine the 
effects of different types of loading. Firstly, it is clear that forces exerted by fixed con- 
straints (pin joints, etc.) are in general included in the theory of Sec. 2, for even though 
such forces may alter in direction and magnitude during a trial displacement, the appro- 
priate component of u,; at this point of application is zero. If secondly, however, some 
of the surface tractions are not of the dead-loading type, then additional terms must in 
general be adjoined to Eq. (15). We consider here only the practically most important 
such forces, viz., pressure-type forces, for which the force applied to a given portion of 
the surface of the body varies in such a manner as to remain always perpendicular to 
that portion and so as always to maintain the same magnitude per unit area. Further, 
the system is still assumed conservative, so that the total work done by these pressure 
forces is independent of the path. If then a pressure P acts on a portion Sp of the surface, 
the work done in the trial displacement u, can be calculated by allowing the intermediate 
displacement to grow at a constant rate—i.e., if ¢ is time, let the displacement at time ¢ 
be (u;t) and calculate the work done from ¢ = 0 tot = 1. This work, w, , is 


W,= [ a { «as» 4 (us) -(—P)-(nd. 


where the subscript ¢ denotes evaluation at time ¢. But 


O(x; + u;t) O(a, + u,l) i: ee 
OX, OX, 





(n,).(dSp), = ; CijKCrng 
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so that 
a ie 7 
W, = | dt | dS; . Cij2lrpal Sip tH Uj pt} { beg Ue at}n, |(—Pu,), 


5 esi 
(25) 


- ] 
— —P | Sp} Witte +, 9 (NU, ku; — MyUy , U;) T 6 Ci KO rpqll; ple, qi ° 


ny ~ 


In this exact expression, the first term would already have been included if P had been 
treated as a dead load; consequently the additional work done is that due to the remaining 
terms. Again we omit terms of third order in u; [see Eq. (10)], and remembering that a 


factor of —2 was incorporated into the derivation of Eq. (10), the term that must be 


adjoined to Eq. (10) is 
| AS pP[njuy uy — Nl, U;]. (26) 


ny 


Considering again the problem of Sec. 6, the term to be added to Eq. (23) is easily seen 


to be 


| (w? + wo’ — vw’ + v’) dé 


and instead of Eq. (24) we obtain 


which is the conventional result. 

8. Adjacent-equilibrium-position method. 
form the condition that an adjacent equilibrium position should exist. Using the notation 
of Sec. 2, and denoting quantities in the perturbed state by primes, the stresses following 


It is now proposed to set up in analytical 


the virtual displacement u,; will be 


: (au \/ Ox! dx? - 
,= p(—) SS, (28) 
\ONng/ OA, OA, 

where x; + u,;. The partial derivatives of U will as before be calculated at constant 
entropy. Expanding the energy term in a power series gives 

‘au \’ ou ou sa l oU - e 

, } a T —. Ors t 9¢ ar ‘ ONrs ONim + a eet 29) 
ONpo/ OMe 9a ONre 2 ON ONre ONim 

where the partial derivatives on the right-hand side of Eq. (29) are evaluated for state IT, 
and where 6n,; is given by Eq. (8). Inserting this result into Eq. (28), replacing xj by 
xz; + u; , using Eq. (3), and neglecting all higher order terms (a process which by virtue 


2 and 3 is considered legitimate) gives eventually 


of Secs. 2 : 
al ; ; aU Ox; OX; OX, OX an 
~ T ij , T li te a T,jU; ee ae 4 ‘ - 7 €, o ’ (30) 

Onn ONr, \OAy OG, OA, Oa, 


where a comma indicates partial differentiation with respect to x; . Then using 
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(within higher order terms) the equation 


0 
_ 5 g M = ( 
Ox’ (ri;) + pl ) 


becomes 


aU ax; dx; dx, Ax 
25 Mies es ———- oo hh = (. (31 
ce FE Tetisrs T [ Onn ON, OA, Oa, Oa, Oa, *' | ) 


In addition to Eqs. (31), a boundary condition must be satisfied. For the case of 
dead loading, the condition is that 


Ti dS’ = T, dS 
and using the fact that 


, Ox, ; 
ni dS’ = (p/p’) ay! n, aS 
Ox 


‘ 
gives as this condition 


( al’ ax, Ox, Ox, Ox . 
: st iy Jn, = (. (32) 


\ On;:; ONps OA, OA; OA; OQ; 
Similarly, the boundary condition for that part of the surface where pressure forces act is 
Ti dS’ = —Pnj dS’, 


whence it is found that Eq. (32) should be altered for this portion of the surface by 
adding to the left-hand side the term 


P(5,,U,,5 — Ue.r)Nq - (33) 


9. Relation between the two methods. It has been remarked in Sec. 1 that the 
methods of Sec. 2 and Sec. 8 can at best be equivalent only for special situations, such 
as at points where an originally stable structure first becomes unstable. Consider therefore 
a structure which follows some stable equilibrium path as the load alters. The path 
will remain stable as long as the second order variation in potential energy is positive 
definite (vanishing only for zero displacement). Consequently, trouble can occur only 
at points where this second order variation [essentially the left-hand side of Eq. (10)] 
vanishes for non-zero displacements. Such a situation of neutral stability will in practice 
be followed by unstable equilibrium states as the load is further increased (see Poincaré 
[17]); we therefore investigate the condition under which the left-hand side of Eq. (10) 
first vanishes for non-zero displacements. Since it always vanishes for zero displacements, 
this condition is equivalent to requiring the minimum of Eq. (10) to be attained at 
non-zero u; as well as at zero u, , and this eigenvalue problem is that obtained by setting 
the first variation of Eq. (10) equal to zero. 

The result of doing this is easily seen to be the same as Eq. (31), with the natural 
boundary condition (32). If pressure forces act on a portion of the surface, the result of 
a variation of Eq. (26) must be adjoined to the natural boundary conditions. This is 
not quite straightforward, because the condition that the pressure loading be conserva- 
tive has not been explicitly stated (without such accessory conditions, the exact pressure 
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work would depend on the path). The appropriate condition turns out to be that on 
the boundary of that portion of the surface where pressure forces act, 


[ ecnst du, dx, = 0 (34) 


for all variations 6u, . This condition would for example certainly be satisfied if the 
pressure surface completely enclosed the body, for then the path length would vanish. 
Alternatively the restraints on u; may ensure satisfaction of Eq. (34)—as in the case of 
a hemispherical shell set on a rigid frictionless plane surface and subjected to external 
pressure. 

Using Stokes’ theorem on Eq. (34) gives 


a 


| ASp[(u,z du,),, — (U, duy) jm = 0 


and if this is used in calculating the variation of Eq. (26), Eq. (33) is indeed obtained. 
It thus follows that at points where Eq. (10) first attains the minimum value of zero 
for non-zero u; , an adjacent equilibrium position exists, and in this sense the two methods 


are equivalent. 
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DISPERSION OF MASS BY MOLECULAR AND TURBULENT DIFFUSION: 
ONE-DIMENSIONAL CASE* 


BY 
B. A. FLEISHMAN** 
Applied Physics Laboratory, The Johns Hopkins University, Silver Spring, Md. 


1. Introduction. If a solute is placed in a solvent in turbulent motion, it is dispersed 
both by molecular and turbulent diffusion. We derive here, in a one-dimensional case, 
a formula for the mean concentration of solute, as a function of z and ¢, in terms of its 
initial distribution, the coefficient of molecular diffusion, and the statistical charac- 
teristics of the turbulent velocity field. 

The one-dimensional problem for the infinite domain is treated as follows. Using the 
initial condition expressed by Eq. (2) below, an initial value problem (for the concentra- 
tion of solute) is formulated for the diffusion-convection differential equation, Eq. (1), 
which contains a rather arbitrary convection velocity v(z, t) which is a function of z 
and ¢. The solution is obtained in the form of a perturbation series, by perturbing with 
respect to the “magnitude” of the convection velocity, and the nth order term of this 
series involves, besides the initial data and the coefficient of diffusion, an n-fold product 
of convection velocity factors, each evaluated at a different point and time. This solution 
is found to be valid at small dispersion times or for small intensities of turbulence. 

Now let the convection velocity be a random (though still continuous and sufficiently 
differentiable) function of x and ¢; that is, let it be a member of an ensemble of functions, 
where the ensemble represents the turbulent velocity field. (This is the approach taken 
in the mathematical theory of turbulence.) When an ensemble average is taken of each 
term of the previously obtained perturbation series, there results a power series repre- 
sentation for the mean concentration in which the nth order term contains the nth order 
correlation coefficient of the turbulent velocity field. This technique, of introducing 
random functions into the theory of partial differential equations, is not a new one. 
Kampé de Fériet, for instance, has treated several classical initial value problems with 
random initial data (see [1, 2, 3]’). 

The initial value problem is solved in Sec. 2. In Sec. 3 a physical interpretation is 
given of a sufficient condition for uniform convergence of the perturbation series solution. 
In Sec. 4 turbulence is introduced, in the manner indicated above, and in Sec. 5 the 
following example is considered: dispersion at small times for initial distributions and 
space correlation coefficients of Gaussian type. 

The investigation described in this paper is now being extended in two directions. 
On the one hand, the treatment presented here will be applied to the corresponding 
problems in two and three dimensions. On the other hand, the formal operations em- 
ployed here, especially those involved in the probability approach to turbulence, must 
be justified. 

The author is extremely grateful to Professor J. Kampé de Fériet for suggesting this 
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line of inquiry. Also, he wishes to thank his colleagues, O. J. Deters and R. J. Rubin, 
for several helpful suggestions. 

2. Solution by perturbation of the initial value problem for the diffusion-convection 
differential equation. The one-dimensional equation governing molecular diffusion in 
the presence of a convection velocity is 


a a's _ _ ale) - 


a 
(see [4], vol. II, p. 593), where s(x, 4) = concentration (mass per unit length) of dis- 
persing matter, D = (constant) coefficient of molecular diffusion, v(z, 2) = convective 
velocity, assumed to depend on z and ¢ in such a way that v, dv/dx and 0’v/dz’ are 
continuous in z and ¢ together. 

For later purposes we assume further that v(z, ¢), v,(z, t) and v,,(z, t) are uniformly 
bounded in an infinite strip -7» <z2<0,0 <¢t< t,/i.e., that there are finite positive 
numbers v”, v> , and v.. , which are the least upper bounds of | v(z, #) |, | v.(z, ¢) | and 
| v,-(z, t) |, respectively, for -~m <r< wo,0St<ch. 

Consider the initial value problem consisting of Eq. (1) and the initial condition 

s(x, 0) = f(x) (-ow< 24 <~o~) (2) 
in the domain —» < z< ©,0 <# < é,. The function f(z) is a given non-negative- 
valued function of z, twice continuously differentiable and uniformly bounded for 
—ao < — < o, We introduce the dimensionless quantities 

r=1T, m= b/T, §=2/(TD), of, 7) = (TD)'*slz, 0/Q, 

(3) 
o(&) = (TD) f(~)/Q, Q2=(T/D)'V, ow, 7) = v(2, d/V, 

where 7’ is some characteristic time for the diffusion process, Q = f=. f(x) dz is the total 

mass present in the system at ¢ = 0, and V is some measure of the magnitude of the 

convective velocity. In the case of a specific (i.e., non-random) convection velocity 

v(x, t) we could let V = v™, while in the turbulent case we shall let V equal the root 


mean square turbulent velocity. 
In terms of these dimensionless quantities, the initial value problem becomes 


Oo 0’o 8(cw) 

— =—- = = C— 1 — —©o@ co a 
- e ge (eS <2, 0S 7S 4) (4) 
of,0)=¢€) (-—e%<§<o). (5) 

It follows from Eq. (3) and the definition of Q that 
[© a =1. 6) 
Assuming that the solution of (4) is expressible as a power series in Q: 

o(€, 7) = 2 MoE, 2), (7) 


we find, after substituting this series for o(é, r) in (4), collecting all terms on the left 
hand side of the equation, and setting the coefficient of each power of 2 equal to zero, that 
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doo _ doo 





ar — ae - 0, (8) 
do, _ Go, _ _ = hes 
Qo = oF = ae (wo;-1) += 1,2, (9) 
The solution to (8) satisfying the initial condition (5) is 
ot, ) = [ o@)TE -&, 2) ae’, (10) 
where 
T'(€, 1) = (4n7)""” exp (—€'/47) (11) 


is the so-called fundamental solution of the heat-conduction equation (8). Furthermore, 
oo(¢, 7) satisfies the normalization condition 


«@ 


/ _oult, 2) de = 1. (12) 


It can be shown, by the use of Duhamel’s theorem (see [6]) and an integration by 
parts, that when w(é, 7) is uniformly bounded in —-» <§< »,0 <7 < 7%, a solution 
of (9) which vanishes at r = 0 is given by 


tr 1 , ~ = 
o,(é, 7) =— [ ee | WO ;-1€ a da, (13) 


where w and o,-_, are evaluated at [ + 2(r — r’)'”a, r’]. Thus when op and o;_,(¢ = 


1, 2, ---) are given by (10) and (13) respectively, the infinite series (7) solves the initial 
value problem consisting of Eqs. (4) and (5) (provided this series and those formed 
by differentiating it term by term converge uniformly in the infinite strip under con- 
sideration; this question will be discussed in the next section). 

If in (13) we introduce the new variable of integration &’ = § + 2(r — r+ 


a )=— fdr [of Mone ETE-e,7- dv, (1h 


\'?a, we get 


where I'(é, 7) is given in (11); (0/0 means differentiation with respect to the first argu- 
ment). Then since [ vanishes at § = +, it can be shown, by inverting the order of 
integration, that f=. o,(, r)d§ = 0 fori = 1, 2, 3, --- . This result, together with Eq. 
(12), implies that solution (7) conserves mass. 

3. Sufficient condition for uniform convergence and its physical interpretation. It can 
be shown that series (7) and the series obtained by differentiating (7) term by term 
converge uniformly in the strip —» < § < »,0 <7 < 7») when 


4(7*) “of liu.b. [| w(t’, 7”) | + | we(€’, 7”) | + | weel€’, 7’) i} re (15) 


—a<t.<@ 
OsT’<S<Te 


Using Eqs. (3) to translate this condition into physical variables, we get 4(é./Dx)'(v™ + 
v™ + v™) < 1, where »”, v2 , and v,2 are defined following Eq. (1). 

This means that the power series is a valid representation of the solution of the 
initial value problem in the time interval 0 < ¢ < ¢ either (1) when for & given arbi- 
trarily, the bound on the convection velocity and its first two z-derivatives is sufficiently 
small, or (2) when for an arbitrary uniformly bounded convection field (twice continuously 
differentiable with respect to x), t is sufficiently small. In connection with case 2 (antici- 
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pating the next section by taking a turbulent velocity field as our convection field) it 
has been noted in [5], p. 99, in considering dispersion from a point source, that “the 
effects of molecular agitation on the dispersion are not always negligible as compared 
with the effects of turbulence; indeed, when the dispersion process starts, the former 
effect is greater than the latter.” 

4. Application to turbulence. Now let w(t, +r) be a random function (though still 
possessing all the regularity properties with respect to ¢ and 7 previously assumed for 
our convection velocity) with (w(t, r)) = 0 at every point (&, 7). The case of a uniform 
non-zero mean velocity, (w), can be reduced to the present case by introducing a new 
space variable, 7 = £ — (w)r. The angle brackets denote an ensemble or stochastic 
average. Thus we regard the convection (i.e., large-scale motion as compared with the 
molecular agitation) as arising from a turbulent velocity field with zero mean velocity. 
It is further assumed that (w’(é, r)) is a constant independent of ¢ and r. For convenience 
we set (w) = 1, which can be done by letting V = (v’)’”’, so that Q is now proportional 
to the root mean square turbulent velocity. (See the discussion following Eq. (3).) 

Rather than o itself, we are now interested in the mean concentration, (c), which 
we propose to find by averaging each of the terms in the series (7). Since oo , given by 
(10), does not contain w, then (oo) = oo . From (14) 


ata=— [def o@, oot’, 21-7 - rae, 
so that c, is linear in w, and since (w) = 0, it follows that 
(oi, 7)) = 0. (16) 
Again, if we use (14) and the preceding expression for o, , we have 
Fi no 


ot, ) = [dr [ we, ETE 8,7 — v1) ae 


x fo de” [ole o)a0le”, 2”) 5g PE — 8%, 1” — 2°) de” 


Thus, oc; involves a product of w’s evaluated at different times and points and there- 
fore, when averaged, will contain the second-order space-time correlation coefficient 
of the turbulent velocity field. Introducing this correlation coefficient: 


R.(é: » T1 5 £2 ’ T2) = ah shes ve _ (w(E, ? 7:)w(E, ’ T2)), 





we get from the expression above for a; , 


(o,(¢, 7)) = [ dr’ = Te =%,1~ 78 i de” 


-o 


(17) 
“i , M7 47 v7 0 , v7 , 4,7 v7 
x | Ra’, 0°58, 2"oolk", 1”) 5g TQ — 8", 7 — 2") de”. 


Similarly, the expression for o, involves a product of n w’s, each with a different 
argument, and so (c,) will involve the nth order correlation coefficient of the turbulent 


velocity field 
R( »71 563,72; °°° 3 &e » Tn) = (w(E, ’ 7:)w(E> ? T2) <2 w(E, ’ Tn) 
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Thus, when the power series representing (c) is known to converge, a knowledge of 
all the correlation coefficients of the turbulent velocity field yields (c) exactly as a function 
of — and r, while if one knows the correlation coefficients R,, up through order n one has 
(in principle, assuming that all the integrations can be carried out) an nth order approxi- 
mation to (c). 

5. Example: Dispersion at small times for initial distributions and space correlation 
coefficients of Gaussian type. We have seen that one situation in which the perturbation 
series (7) actually solves the initial value problem is when, for a given convection field 
(or class of convection fields), the time interval 0 < 7 < 7» is sufficiently small. Let us, 
then, calculate expressions for a(t, 7) and (o2(&, 7)) at small values of 7, for certain 
initial distributions and correlation functions. (We recall that the first-order term,(o;), 
vanishes.) In terms of the computation the restriction to small values of 7 is imposed 
by the desire to approximate certain functions of 7 so that they can be integrated (see 
Eq. (22) below). 

Consider initial distributions of the form 


o(t) = (4er*)'” exp (—£°/4r*), 


depending on the parameter r*. Then from (10) 


oof, 7) = [4e(r + 7*)]-'” exp [—€/4(7 + 7*)]. (18) 
(This is the distribution that would result from molecular diffusion of a unit mass con- 
centrated at ¢ = 0, starting at time r = —7*.) 


We assume that the turbulence is homogeneous and stationary, and furthermore 
that the second-order correlation coefficient can be written as a function of the ¢’s times 


a function of the 7’s: 
Rae’, 1758", 17) = Oa — 85 7 — 7") = AE — EB — 7”). 


Finally we consider the particular class of functions 
R-() = exp (—£°/4), (19) 


with the parameter b. Without specifying ®, , we assume that in the neighborhood of 

= 0 its rate of decrease is sufficiently slow so that for the small values of r under 
consideration &,(r) - 1 (since ®,(0) = 1). We shall, however, carry along the factor 
®, in the calculations until it is necessary to invoke this assumption, so that perhaps 
some enterprising reader can find an explicit functional form for ®, which will obviate 
the approximations which are made at that point. 

Under the foregoing assumptions, if we replace ¢’’ and r’’ by the new variables of 


integration §, = : — §’’ and rz = r’ — 1’’, Eq. (17) can be written 





(.(, r)) = — 4? mf es a G—ry" ‘. (¢’ — &) ex | - $5, oe a dg’ 
(20) 
x val ; . Seley ors [ ERe(E2) oo(E’ —&,T — 72) exp (- J g: =) dé, , 


where 9I'/d¢ has been obtained from (11). Inserting the expressions for a, and ®; given 
by (18) and (19) into the first integral to be evaluated, 


Ii’, ',%) = [. ERe(Es)oolt’ — & , 7’ — 7) exp (- &) dg, , 
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we find, after an essentially simple integration, that 


(2 LOE op [ln tO] 
lj _ D®”? exp 4D . 


where 
D = 122’ — 12 + 1*) + Wr’ + 7°). (21) 


Next, we must evaluate 


dn “ee 
Pe sm [ R,(72) E(t. + »] 29 
~ 4r J, D™ exp | 4D drs , (22) 


where D is given by (21). Since r, < r’ < r <1, R,(r2) & 1. Then, by considering 
different relative magnitudes for 7, r* and b, we can simplify D and evaluate J" in three 
cases: (1) r <b, r* Kb; (2) r Kb, r K 1*; (3) tr K r*, Db K 1*. In each case we shall 


evaluate J’ (é’,.r’), then (see Eq. (20)) 


ts» r\ 2 , (¢’ — §)’ Ther , , ‘ 
He) =f @ -dew|- ya Ire, oe, (23) 
Was 4(r — 7) 
and finally 
— a ey i(é, 7) +, 2 
(o,f, 7)) = — 4a” | io = ry" dr’. (24) 


Case 1:7 < b, r* <« b. From (21), neglecting terms of second order in the small 
quantities 7, 7’, 7, and r*, we get D = b(r + 7*). We insert this expression for D and 
R, = 1 in (22) to obtain 


"i dr, 


, > 8?) edn 
J" = Ar 2 i exp | Ar’ .s =| (_ > 7)? 





—— a se | 
= 8! + ry? OP | a +) I 


Then, using this result in (23) and integrating gives 


7 7! f Wil a (’ — §)° as | ie! 
; = 7 (@’ —) eaxp| - er -—- Sr er I 
ie (4 + rt)’ Jew " 4(r — 7) 4(7r’ + 7r*) J” 


re? — 2(r + 37*) P : 
(o.(€, 7)) = & aa at) exp | - | (25) 


=" 16n°”° (7 + r*)”/ 4(r + 1*) 
Case 2:7 Kb, r K r*. Now D = br* and 
r ¢/ at’ \ 
v= — os | exp (- 3) dr2. 
4g sg? Jo ' ér° c 
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Proceeding as in Case 1, we find 


, Fe Dig 2 
(ult, 9) = gga an exp (— 5) (26) 


Case 3:7 K< 7r*, b Kr*. Now D = 1*(r. + b). By performing the successive integra- 
tions indicated above (and using along the way the fact that r’ << r < 1r*), we obtain 
in this case 

_ * 2 
(o2(é, T)) = ane exp (.. £ on singe 20° [(r + b)'”? neal oy} (27) 

Remarks. 1. The expressions for (c,) in Cases 1 and 2, given by Eqs. (25) and (26) 
respectively, are quite similar in form; indeed, if in Case 1 we made the more restrictive 
assumption that r < r* < b, (25) would reduce to (26). Equation (27), the result for 
Case 3, however, is different from the other two: only in this case does b appear in the 
expression for (¢,). Thus (Az) = oo + 2*(e,), the second-order approximation to (¢2), 
does not depend on b when r < b, whereas in the absence of this assumption the expres- 
sion for (A,) does involve b. (It should be noted that b is proportional to L’, where L 
is the length scale of turbulence, the integral of the correlation coefficient from 0 to © 
with respect to 2.) 

2. If in Case 3 the more restrictive assumption that r << b < r* is made, we get 
(o,) = 0, for 


1/2 
f 1/2 1/2 /2 /2 ef 
(r + b) — >’? = p' (1 + r) = 1] ~ b' (1 ~ =) = | - ap * 


3. A comparison of Eqs. (18) and (25) shows that in Case 2 (¢,) = (r°/2)d0o/dr. 
It is not clear what significance, if any, there is in this interesting relation. 

4. Within the limitations of requirement (15), our sufficient condition for the validity 
of the perturbation scheme, it can be shown that in the region of the &, 7-plane where 


0.06-— 
Approaches €-axis 


asymptotically 





0.04} —_ 
| | | | | 
| | 


0.02} 


0 } 


<%>/T* 


-0.04| 


-0.06| 





“0.08 5-4 2 - oO | 2. 3 5 6 





6 


Fic. 1. Example of second-order correction to mean concentration at small times: (o2)/r? plotted as a 
function of ¢ from Eq. (26) with r* = 1. 
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(A;) is a good approximation to (c), 2’(c,) is much smaller than a, . (This is not surprising 
since in perturbing with respect to the convection term, we assumed this term is small 
in comparison with the others in the original differential equation.) It should be kept in 
mind, of course, that sufficient conditions for a certain result flow out of the particular 
mathematical technique employed, and often the result is valid under much broader 
conditions. This is particularly true of perturbation methods. Nevertheless, even if we 
limit ourselves to values of the parameters and variables which satisfy (15), the results 
are of interest. Equations (25), (26) and (27) show that in all three cases, for given 
values of b and r* and a fixed value of 7, (c,) is negative at ¢’ = 0, increases, with in- 
creasing £’, to a positive maximum, and then approaches 0 asymptotically as #? >. 
This is illustrated in the figure for Case 2. The exponential factors in both o, and (a) 
make them go to 0 as ” ~~, but the relative magnitude of (o,), as compared with ap , 
increases with £; indeed, (¢,)/¢, becomes proportional to ~’ as ~ increases. Adding 
turbulence to molecular diffusion, then, ‘‘sweeps out” the initial distribution more 
rapidly, and this additional dispersive effect assumes a particular functional form, at 
least in this set of examples. 
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WAVE PROPAGATION IN RODS OF VOIGT MATERIAL AND VISCO-ELASTIC 
MATERIALS WITH THREE-PARAMETER MODELS* 


BY 
J. A. MORRISON 
Brown University 


Abstract. The problem of an impulsively applied and subsequently maintained 
constant velocity at the end of a semi-infinite rod of Voigt material is considered and 
integral expressions are obtained for the velocity and stress distributions in the rod. 
Apart from a constant factor the stress distribution arising from an impulsively applied 
constant stress at the end of the rod is the same as the velocity solution above. The 
same problem is considered also for materials with three-parameter models. The stress 
distributions for both the case of constant applied stress and constant applied velocity 
are represented graphically for the materials considered, dimensionless co-ordinates 
being used. 

1. Fundamental considerations and the Voigt model. Lee and Kanter’ considered 
the problem of an impulsively applied constant velocity on the end of a rod of Maxwell 
material and analytically determined the subsequent stress distribution in the rod. 
Glauz and Lee’ treated the same problem for a four-parameter model visco-elastic 
material and determined the stress and velocity distributions by the method of charac- 
teristics. In order to complete the picture we consider here the other two-parameter 
model, which corresponds to the Voigt material, and also the three-parameter models. 

The Voigt model is shown in Fig. 1 and is composed of a spring and a dashpot in 

E 
‘sume A'4 4 44 foe | 


€ 

















vel | 
Fia. 1. Model for Voigt material. 

parallel and so the Voigt material presents a retarded elastic response to stress. The 

stress strain relation is 


cus, (1) 
m 


where o is the stress, ¢ the strain, the suffix denotes differentiation with respect to time 
and E and y are elastic and viscous constants, respectively, for the material. 

Let us consider a creep test, together with unloading, in which a constant stress oo 
is suddenly applied at time t = 0 and suddenly removed at ¢t = t, so that o’ = oo{H(t) — 
H(t — t)} where H(t) is Heavisides’ step function, 

1,t>0 
H(t) = 
0,t <0. 
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From Eq. (1) it is found by integration that « = {y()H() — y(t — t.)H(t — t)}, where 
¥(t) = (o)/E)(1 — e*“'). For convenience we introduce dimensionless variables r = Eut, 
>’ = o'/o, and & = (Ee)/c, . Then, 


2’ = {H(1) — H(r — »)}, 


5 = { o( 7) H( tT) — aT = to )H(r = To)}, 
where 
g(t) = (1 — ee’). (3) 


The creep test is represented in Fig. 2 where r> = Et, has been given the value 2. 





We now turn our attention to the problem at hand. We are concerned with the propa- 
gation of longitudinal waves in a semi-infinite rod, x > 0, which is initially unstrained 
and at rest. Let u(x, ¢) denote the displacement of the sectton xz of the rod so that its 
position at time ¢ is given by (x + u). Then, if p is the density of the unstrained material, 
(4) 


(5) 


Pee mes 5 

«€e= —@, ; 

where suffixes denote partial differentiation with respect to the corresponding variable. 
Here o and « are the nominal compressive stress and nominal compressive strain, respec- 
tively and it is supposed that the stress strain relation of Eq. (1) applies to these nominal 


values. 
Eliminating o and ¢ from Eqs. (1), (4) and (5) we obtain 


; l : 
pure = Euz, + i Uret ] (6) 
be 


from which it is readily verified that the stress o and the particle velocity v = u, both 
satisfy the partial differential equation 


hs l “ 
Pf = Efsz + — fast - (7) 
Le 


It follows that the stress distribution o’(z, ¢) in the case in which a constant stress a 
is suddenly applied at the end of the rod and then maintained is of the same form as 
the velocity distribution v(z, ¢) for the case of a suddenly applied and subsequently 
maintained velocity vp at the end of the rod. In fact, 


v(x, t) o’(x, t) 
wh Seid a eee (8) 
Vo To 
If we introduce the dimensionless variables § = (pE)'/*ux, V = v(z, t)/vp and, as before, 


r = Ent, >’ = o'(z, t)/oo , Eq. (8) becomes 


z'(é, r) = VéE, 7). (9) 
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We therefore confine our attention to the problem in which the end of the rod (which 
is initially unstrained and at rest) is suddenly given a velocity vp which is subsequently 
maintained and determine the stress and velocity solutions, o(z, t) and v(x, t). We 
introduce the dimensionless stress 2(¢, r) = a(x, t)/(pE)'’’v, . We will then have deter- 
mined the stress solutions Z(t, 7) and 2’(t, 7) both for the case of constant applied 
velocity and constant applied stress. 

The determination of the solutions, for the Voigt material, is carried out in Appendix 
A by the method of Laplace transforms. Zverev’® adopted the same method but carried 
out the inversion of the Laplace transform in a different manner and obtained an infinite 
integral expression for the velocity distribution. The present author gave a similar 
expression for the stress distribution‘ but the method adopted in this paper is readily 
extended to the case of the three-parameter model with two viscous elements and in 
addition the finite integral expressions obtained are readily computed for the moderate 
values of 7 in which we are interested. 

The stress solutions £(£, r) and D’(£, r) as given by Eqs. (A.18) and (A.14), respec- 


~ 


tively, are depicted graphically in Figs. 3(a) and 3(b). In Fig. 3(a), = is plotted against 
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Fic. 3(a). Voigt model. 


t for different values of r so that the distribution of stress in the rod at various times, 
arising from an impulsively applied and subsequently maintained constant velocity 
at the end of the rod, is illustrated. From Eq. (A.19) it is seen that the stress at the end 
of the rod, which is very large immediately after impact, decreases rapidly at first and 
subsequently continues to decrease, but with diminishing rapidity, ultimately approach- 
ing the value unity. The stress at a position of the rod very close to the end rises extremely 

3], N. Zverev, Prik. Mat. Mek. 15, 295-302 (1950), (Russian). (Brown University translation 
A11-T12/16). 

‘J. A. Morrison, Wave propagation in a rod of Voigt material, Tech. Rept. PA-5/17, Div. Appl. Math., 
Brown University, Providence, R. I. (1953). 
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rapidly shortly after impact to a large value and then commences to decrease. The 
occurrence of large stresses near the end of the rod is due to the response of the viscous 
element to an attempted sudden change in strain. Such a change in strain would cor- 
respond to an infinite strain rate and so demand an infinite stress to produce it. 

However, as the time increases the viscous behavior tends to die out and the rod 
tends to behave as in the case of a completely elastic material. At each position in the 
rod the stress ultimately approaches the value unity. It is noticed that there is a con- 
siderable spread in the stress distribution, as opposed to the sharp wave-front occurring 
in the case of an elastic material, and that the effect of the disturbance is propagated 
instantaneously (in a decaying manner) to all positions in the rod. 


r’ 
1.0 
9+ 











Fic. 3(b). Voigt model. 


In Fig. 3(b), 2’ is plotted against ¢ for different values of 7 so that the distribution of 
stress in the rod at various times, arising from an impulsively applied and subsequently 
maintained constant stress at the end of the rod, is illustrated. At each position in the 
rod 2’ ultimately approaches the value unity but it is again interesting to note the 
considerable spread in the stress distribution. 

2. The three-parameter model with one elastic and two viscous elements. There are 
four true three-parameter models, all the others reducing effectively to models with 
fewer parameters. We first consider the two models shown in Figs. 4(a) and 4(b) which 
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Fia. 4(a). Fia. 4(b). 
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have two dashpots and one spring. For the model of Fig. 4(a) we have 








Ee’ + Ve =g¢g= 2G. — é). (10) 
mn mn 
Hence, eliminating ¢’, 
Ee, + “+ = Evo + (1 + te : (11) 
mn B 
In a similar manner we obtain the stress strain relation for the model of Fig. 4(b), 
* 
Pa + p*o = a(t + us) + ron 5 (12) 


If we set (E/E*)'? = p*/u = u*’/u’ = (1 + w/u’), then Eqs. (11) and (12) become 
identical. Consequently we confine our attention to the model of Fig. 4(a). If we set 
o = ooH(t) in Eq. (11) and integrate we find that 


e = 2 (But + (1 — HO, (13) 
so that the creep test with unloading is given by Eq. (2) with 


or) = {r+ (1—e™”)}, (14) 


where A = u/p’. This is depicted graphically in Fig. 5 for the particular case \ = 1. 








Fra. 5. 


It is noted that after the removal of the applied stress the strain approaches an asymp- 
totic value given by « = 7») (= 2) and this residual strain is due to the viscous flow in 
the dashpot with viscous constant » which is in series with the retarded elastic or Voigt 
element. The strain of the delayed elastic portion is asymptotically recovered, as is 
seen from Fig. 2. 

Eliminating o and ¢ from Eqs. (4), (5) and (11) and setting v = u, , 


pE pv, + o(1 + raat _ Ev,, + ar Yas ’ (15) 


and it is readily verified that o satisfies the same equation. Hence Eq. (9) again holds 
and the stress distributions =(é, r) and =’(€, r) arising when the end of a rod, initially 
unstrained and at rest, of a material corresponding to the three-parameter model of 
Fig. 4(a) is impulsively subjected to constant applied velocity and constant applied 
stress, respectively, are determined in Appendix B, Eqs. (B.8) and (B.14). 
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~ and >’ are plotted in Figs. 6(a) and 6(b), respectively, against ¢ for different values 
of 7 in the particular case \ = 1, so that the two viscous constants in the model are equal. 
Thus in Fig. 6(a) the stress distribution in the rod at various times, arising from an 
impulsively applied and subsequently maintained constant velocity at the end of the 
rod, is illustrated. The stress at the end of the rod is very large immediately after impact 
due to the combined ‘“‘rigid’’ responses of the viscous element and delayed elastic element, 
in series, which comprise the model of Fig. 4(a). Nevertheless it decreases very rapidly 
and has a value less than unity at time r = 1. The stress at the end of the rod continues 


z 
1.0} 














Fic. 6(a). Three parameter (1 elastic, 2 viscous). 











Fic. 6(b). Three parameter (1 elastic, 2 viscous). 
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to decrease as 7 increases and at each position in the rod the stress ultimately vanishes. 
In Fig. 6(b) the stress distribution in the rod at various times, arising from an im- 
pulsively applied and subsequently maintained constant stress at the end of the rod, is 
illustrated. At each position in the rod 2’ ultimately approaches the value unity. 
3. The three-parameter model with one viscous and two elastic elements. We now 
consider the remaining two true three-parameter models which are shown in Figs. 7(a) 
and 7(b), these having two springs and one dashpot. 























(o-0" E (o-o 5° 
é spiel ita 
"af 
oT . o 6 : o 
, rT 
Cg g | L1/pa* 
rss 
Fic. 7(a). Fic. 7(b). 
For the model of Fig. 7(a), 
, 
é = a + po’, 
(16) 
l ( " 
=s — o'). 
€ E og 
Hence, eliminating o’, 
t E " Ps 
_ + yuo = e(1 + z) + Eue. (17) 
In a similar manner we obtain the stress strain relation for the model of Fig. 7(b), 
‘ 1 E* CO; 
E*e + —*™ o(1 + E) + BY a (18) 


If we let E*/E = E*'/E’ = (u/u*)'” = 1 + (E£/E’), then Eqs. (17) and (18) become 
identical so that we confine our attention to the model of Fig. 7(a). If we set ¢ = ooH(t) 
in Eq. (17) and integrate we find that 


oo E’ EE'pt } 
=F a OF — ee ae , 19 
€ E {1 (E+ BE’ exp | (E +E’) H(t) (19) 
so that the creep test with unloading is give by Eq. (2) with 
1 —T o 
= <1 — ———- ep | -——_ = 20 
= ft a+ la + al ane 


where k = E/E’. This is depicted graphically in Fig. 8 for the particular case k = 1. 
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For this model there is an instantaneous elastic strain of magnitude o,/(E + E’) = o./E* 
when the stress is applied and the strain instantaneously falls by this same amount 
when the applied stress is removed, and asymptotically tends to zero. (This behavior 
is perhaps more apparent from Fig. 7(b), in which there is a spring of elastic constant 
E*’ in series with a retarded elastic element.) 

Eliminating o and ¢ from Eqs. (4), (5) and (17), 


tte E 
2 + puu., = len I + z) + Eyu,, , (21) 


and it may be verified that o and v satisfy the same equation. The stress distributions 
Z(é, 7) and 2’(é, r) for the problems previously considered for the other materials are 
determined for this material in Appendix C, Eqs. (C.20) and (C.19), respectively. From 
the prescribed condition on the end of the rod 2’(0, r) = 1 and a simpler expression for 
=(0, r) is given in Eq. (C.25), so that these provide a check on the accuracy of the 
computation of 2’(é, r) and 2(é, 7) from Eqs. (C.19) and (C.20). 

> and 2’ are plotted in Figs. 9(a) and 9(b), respectively, against ¢ for different values 
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Fic. 9(a). Three parameter (2 elastic, 1 viscous). 


of 7 in the particular case k = 1, so that the two elastic constants in the model are equal. 
Thus in Fig. 9(a) the stress distribution in the rod at various times, arising from an 
impulsively applied and subsequently maintained constant velocity at the end of the 
rod, is illustrated. A wave of stress discontinuity is propagated along the rod with con- 
stant velocity (2E/p)'” (or in the general case [(1 + k)E/kp]'” = [((E + E’)/e]"”), 
relative to the unstrained rod. Hence at time r the portion of the rod for which > (2)'’r 
(or > (k + 1)'”7/k’”, in the general case) is unstressed. The magnitude of the stress 
discontinuity at = (2)'7, (or & = (k + 1)'77r/k'”), is (2)'” exp (—7/4), (or [(kK + 
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1)/k]'” exp [—1r/2k(k + 1)]), so that it decreases exponentially as the wavefront ad- 
vances along the rod. In the case considered, the stress at the end of the rod jumps to 
the value (2)'”” immediately after impact and then monotonically decreases, ultimately 
approaching the value unity. The stress at each position in the rod also ultimately 


approaches the value unity. 
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Fic. 9(b). Three parameter (2 elastic, 1 viscous). 


In Fig. 9(b) the stress distribution in the rod at various times, due to an impulsively 
applied and subsequently maintained constant stress at the end of the rod, is illustrated. 
Here again a wave of stress discontinuity is propagated along the rod, the constant 
velocity of propagation being the same as before. In this case the magnitude of the stress 
discontinuity at & = (2)'”’7, (or at & = (k + 1)'”7/k'” in the general case), is exp (—7/4), 
(or exp [—7r/2k(k + 1]), so that again it deereases exponentially as the wavefront ad- 
vances along the rod. The stress at each position in the rod ultimately approaches 
the value unity. 

4. Concluding remarks. We have been concerned with the stress distributions in 
rods, of Voigt material and of visco-elastic materials with three-parameter models, 
impulsively subjected to constant applied velocity or constant applied stress at the end 
of the rod. We have seen that there are two basic types of true three-parameter models, 
namely those with one spring and two dashpots and those with one dashpot and two 
springs. In a subsequent paper it is intended to compare the stress distributions =(é, r) 
and >’(é, r) for seven different visco-elastic materials; elastic, viscous, Maxwell, Voigt, 
the two basic three-parameter models considered here and also the four-parameter 
model considered by Glauz and Lee’. With this in mind we chose the model of Fig. 4(a), 
in preference to that of Fig. 4(b), since the ultimate behavior may then be compared 
with that of a dashpot with viscous constant yu. Similarly we chose the model of Fig. 
7(a) since its ultimate behavior may be compared with that of the Voigt model considered 


in Sec. 1. 
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integrations. 
Appendix A. From Eq. (7) it is found that V(é, r) satisfies the equation 
¥ ts = Ver + Veer ° (A.1) 
Also, from Eq. (4), 
Re = Za ° (A.2) 


V(o, r) = 0 = Yo, 7), 
(A.3) 
VO, 7) = 1 
Applying the Laplace transform, 
Lif(é, a} = | e  “F(E, 7) dr, (A.4) 


70 


to Eqs. (A.1-A.3) and remembering that the rod is initially unstressed, unstrained and 
at rest, we obtain 

GLiV} = (14+ QIL{ V3 Jee = (1 + QL{ Ve} (A.5) 

gL} V} = —(Liz le = —L{z, (A.6) 


and 


LiVlo,7)} =O = Lix(o, 7}, 


= l 
Liv, r)} =-- 

q 
Solving the ordinary differential Equation (A.5) for L{ V(é, r)} and applying the boundary 


conditions (A.7), we obtain 


ss 1 —§q | ’ 
} r)} = —~ eX las > ° f 
L{ V¢é, 7)} _ F is (A.8) 
From Eq. (A.6) and the boundary condition in (A.7), 
; (1 + q)'”” —tq | : 
{p(e, 7)} = ——*—_ exp | =“ |. A.§ 
L{ X(é, 7)} ; exp | 73 ai (A.9) 


In carrying out the inversions of L{V} and L{2} we make use of the followin 
; j j g 
general result for Laplace transforms which is valid under certain conditions.” 
If 
F(q) = Li{f(n} and —— exp [— nh(g)] = Lig(r, n)} 
g 


\q) 





6H. S. Carslaw and J. C. Jaeger, Operational methods in applied mathematics, Oxford University Press, 
London, 1947, p. 259. 
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then 
1 7 
dé F{h(q)] = 14 | e(r, n)f(n) an}. (A.10) 
Now, 
1 1/2 xp (—#?/47 
(yi oP [—t(q) “"] = ere i (A.i1) 
and 
hs -" n) a cos [2(n7)'”* i 
(q) 1/2 exp ( q coe Lye, (wr)? 
so that 








to 


= o 1_ 4) ] _ feos [2(n)"r = 07) one \ 
(q)” exp | Aq + = 2) = 14 (nm) — 9)? H(r — n) (A.12) 


Hence, from the result of (A.10), for qg > 1, 


~~ t(q — | _ {! ” eos [2(n)' (7 — n)' - (2 a | F } (A 
(q — 5 [= @g? J~* Ji wes rags S 





Applying the shift theorem we obtain from Eq. (A.8), 


(E , e" f" cos [2(m)'"(r — n)'") (2 ) 
> - T = f 7 =. = “<:. ees —_ & _= —- A. 
>’(é, 7) V(é, 7) > [ Le “ ni? exp |2n i dyn. (A.14) 
Further, 
1/2 
exp (. 2) = ae tT) — (2) J ,(2(n7)' , (A.15) 
q { T ) 


where 6(7) is the delta function. Thus, 


] PP (n)'” F 1/2 1/2 
exp | -a(« + . 2) = Lye E ~ = oon 73 J1}2(n) “(7 — 9) 
- A(r — »} (A.16) 


From Eqs. (A.10) and (A.11) it then follows that, for g > 1, 
(q)'”? —i(q — 1) exp (27 — £°/4r) 
r 7) ©Xp ; ia—|=L (x7)? 


(q — 1) (q) 
J[2)'(r = 0?) | (2 ~f) } nial 
[ (x)'(r om n)? exp |< ’ dn ° (A.17) 











6R. V. Churchill, Modern operational mathematics in engineering, McGraw-Hill Book Co., Inc., New 


York, 1944, p. 299. 
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Applying the shift theorem we obtain from Eq. (A.9), 


_ _e_ Jexp (—’/47) 
z(E, r) = fin (7)? 


7 Js[2(m)'*(r — 0)"*)] exp | -2%s —n- | an}. (A.18) 











J (r — »)'” 4n 
Also, since 
(1 + q)'” (q+1)” 
L{s(0, = = 
120, 7) q (@+i—1)’ 
we find that 
>(0, 7) = cam + erf co}, (A.19) 
and 
@ se a —e 
dr {z(0, t)} ™ (mn) (1)? . 


Appendix B. Equation (15) in dimensionless form becomes 
Ve t(L+A)V., = Vee t+AVee- « (B.1) 
Equations (A.2) and (A.3) still hold and applying the Laplace transform we obtain 
g{l + (1 + A)gIL{V} = (1 + AQLiVex}, (B.2) 


together with Eqs. (A.6) and (A.7). Solving Eq. (B.2) for L{V(é, r)} and applying the 
boundary conditions of (A.7), 











: _1. fs ggfii ++ na” | 
Liv (é, 7)} = exp { (1 rw dg)” . (B.3) 
From Eq. (A.6) and the boundary condition in (A.7), 
(1+ Ag)” { tg{l + (1+ yal 
{> - 5 exp 4 — 72 - : 
L{z€, t)} (q)*”*[1 + (i + rA)q]'” exp (1 + rq)” (B 4) 


The inversion of L{ =(é, r)} will be carried out in a manner similar to that of Appen- 
dix A. 





1 (2+2)\] _ 
-” [x2 + 0+ Ng M+ oh dite 
where 
(7,7) = {exp [e+ MT a — n) 


(n)™* 2(n)'/"(r — a" > } , 
— xd + ») (7 _ n 14 Md + Vis H(r n) |?- (B.5) 
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Hence, using the shift theorem, 


mrgil + (1 + Adg] = -r/h 
0 {— Ha tage oH He 





But, 





(q 
Hence, from the result of (A.10), we obtain 


. a 1 r M1 +r)? E (1+ ne] 
26) = a? LT oh (yt PL aXe 








1 t(1 +f ate he = 1 (1 + NH? 
aie [- 0)? ~ Ney “ - 4yz } 








+ 2m) “(7 — 9)? (2+rA(r—) (+ ne] dn } 
- | rf A(1 + A)” | sa | - A(1 + A) 7 4Xn (r — )'”*)° 


We now carry out the inversion of L{V(é, r)}. 


1 1 (2+A)l] _ 
1 ep [ —a{a + qty - SS} = L{x(r, )}, 





where 








_ n(2 +d) | | 2(m)'(r — 2)” _ 
x(r, n) = exp | 22 + » Jf AI + _ lac 0). 


Hence, using the shift theorem, 


A —nqil+(1+A)q]l _ py -a 
(i + dq) exp { (1 + 0 + 9) } = Lie"'x(7, )}- 








From (B.6) and (B.10) we find that 








N71 + (1 + dg] {=ahalt + (1+ wal} ee ; 
(1 + Aq) - (1 +0 + AQ) Lie"""[A(L + A) (7, 9) 


ies x(r, n)}} ° 


But, 





1 —¢(1 + ro | Ls { [ + wr} 
q exp | (A) 1/2 = L erfc 2(r)'? ° 


Hence, from (B.11) and (B.12), using the result of (A.10), 


a i T, 1+ )'” 
rier = veo ee [oem — 2 | ete [ELE | a, 


(B.6) 


(B.7) 


(B.8) 


(B.9) 


(B.10) 


(B.11) 


(B.12) 


(B.13) 
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Thus, finally, 


-_ ee ee, | +d) 
>’(é, 7) = exp E + i et | Xan? | 
ar f = / 1/2 1/2 92 1/2 -. 1/2 
1p {yfaimte 9] , ae mrig , [xem — 9") 
A1 +A) Jo | wi +" (7 — ») (1 + A)’ J 


~ 2 + d)(r — 2) co | 
exp | =! x1 +d | exte ! an)? dn}. (B.14) 


Appendix C. Equation (21) for the velocity becomes, in dimensionless form, 








Vie tkhViere = Veet (Ll t+ hyve - (C.1) 

Equations (A.2) and (A.3) again hold and applying the Laplace transform we obtain 

g(ii+tk@™LiV} = [1+ 4+ Ag]L{Ve}, (C.2) 

together with Eqs. (A.6) and (A.7). Solving Eq. (C.2) for L{V(é, r)} and applying the 
boundary conditions of (A.7), 


? a 1 ; [ —Eq( 1 + kq) om 
V(é, 7)} at fi + a + Ba)” if (C.3) 





From Eq. (A.6) and the boundary condition in (A.7), 





on i+(d+hq gat + eT 
ue. 07 ~ BESS eo i7a- C. 
l ‘ q(1 + kq)” = ep 1 + k)q]'” (C.4) 
Let 
; {_—tg(1 + kg)? rr r 
exp li + a+ bal ~ L{A(é, 7)}, (C.5) 


so that A(é, 7) is the dimensionless acceleration. Then, 


oT 


VE, 7) = | A(é, ¢) df, (C.6) 
and 
3, 7) = [ A(t, 7) di (C.7) 
Let 
ian 1a 
bl at : = aE FD aaa ee wand 
Then, 





as pate+ Or |. p{- Egil + ka co 
ep | ~« Vqta-p) “PV i++ hq” ” 
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Now, 


exp | -#(e ~~ t # 2 | = exp (—wq) exp [—w(a — 28)] exp | - Ala - 8 


(¢g — B)' 
= _ 28(a — 8) 7” . 
q + B] [(q° am > pany q + 8) . (C .10) 


~ 








Also, 
fexp [—n(q’ — 8°)'”?] — exp (—ng)} 

= 7 Aa 7) 1,[8(7° — 9°)? )H(r — n}. (C.11) 
so that 


apg is Qn)"?}H(r) }>. (C12) 
+ rapes 9 Dn)” I, {B(r)'”? “(r+ 2n)' }  } (C 


From (C.12) and (A.10) it follows that, if F(q) = L{f(r)}, 














F[(q° _ 87)' — q 4. B] = 1A (8) 6(7) 
on, LalB(2)'"(r + 20") \ 
+8fe Oe art Ic) dap. (C.13) 
But, 
1 
exp | - ata _ kk - ah | = 1A a+ — w) 
— (w)'/? a 2(w) foe | ba } © 
+ E+ DG — 0) 1 (e+ 1) JP — oy. 14) 
Thus, since 
28(a — B) = — 
. kk + 1°’ 
[ __28(a — 8) 
ep | - olla — By a4 8 er rat | 
; -sw 1,[B(7) + 2w)'”?] 
- un [w(2e — 38)]8(7) + Bue Se 7 te 





° py 1,[B(2)'"(r + 20)'*] 
sd Oe ae 


(co)*/? 2(es)'/?(9 1s ad dn } , 
te lars la ine 
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Hence, from (C.10) and (C.15), 


exp E —w(qg — a) ( at + a 4 = I\exp [w(a — B)]i(r — w) 


+ Bw exp [— da - 9) =s7 ie mia 





(7 — ww)” 
” BA (r —— a 2)1/2 
+ B exp [—w(a — 6)JH(r — w) ; e F"(n + w) ae ta pte A 
— ue asl dn ats} , 
Cateae outed oie wom 


Thus finally, using the shift theorem, from (C.5), (C.8) and (C.9) 
—  f-2k +17 /f__ ore | | _ _(e% | 
AG, = om | Seep he la + D7 )'L7~ &+" 
a poh | | _ porte 
+ Xk) (hk + 1)” ex | (k + 1? H| + (k + 1)? 
J {7 xml ye 
| é | _@+D 
{eo ye ‘Lo 2k + 1D) 
Ll" (k+ 0 
— —7 (n+ 8G)" i 2(én)'”? | 
+ I “_ Ernee + 7 | (k + 1)°*(n)'” ; (k + 1)” 


{| +} (k)?9 ii _ ke + 
(k + 1)” (k + 1) 


2k(k + 1) 





























7? 


. Ss dn 
{|= ‘ (k)'/?n | S k(é + 7 
{ (k + 1)” (k + 1) 


|* 
»— en | ret De _ ws — _ (eM | 
tT) — ep | Dkk + 1) +)” ate &+ 7) (C.18) 


then from (C.6), (C.7) and (C.17), 


(C.17) 








If we let 
Be, 7) = Aé€é, 





( 
>’, 7) _ Ve, T) _ ‘exp | 37 rT a 1)°7 2 | 


tT (k)'¢ ] . 
B dt? + — =——05 | C.1 
t I g, p sp E (k + 1) ' c 9) 


and 


ee S(k- oe ia x laa | 
Zé, 1; = | (k) \172 €x p Qk(k + 1) 


[(k+1)2/97}/(k)*/9 (k)¢ ] 
+ I B(t, 7) ar\a| + -EGD@} 20) 
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From (C.17) and (C.18) we obtain 








_ (4k +E | —t ] 
oo. ax {BE, 7)} = 8(k)*(k a 1)” exp 2k) *(k 7 1)” 
o>’ 
is 
sigabiniaammmaiaiaiie Or @ ”) ( ) 
and 
: _ (4k +1)r [ —T | 
teeter ee ayy BE 2) = BE + IP LIKE + D 


. dz ‘ 
= lim {- 8 (é, oh. (C.22) 


[E+ (41) 2/97) / (ke) 78] 


Also, from (C.19) and (C.20), 


‘ rp} 
lim ‘2 (é, } 
[rol (e)*/9E)/[ (e4+1)2/9) +) or 














— [4k + DE - 4'7%K + 7) | -t | 
— 8(k)*(k + 1)’ exp 2(k)'"(k + 1 a , (C.23) 
and 
az’ _ [(4k + I)r + 4k(k + 1) | —1 | 
canal yas , of =" 3@7e4+ D7 PLaxE+ DS 
Further, 
. ys = nes ae 
a... (2"G, 9} E= + mi 
(C.24) 
_ (k+1)” p| —t | 
or. ee (2, 7} (k)'” _ 2kh(k + 1) J 


From’Eq. (C.4), 


m _fl_+ad+hq)” 
L{z(0, t)} _ q(1 + kg)” 


1 
[a+n+3] 1 


= (7 > ky? { ‘ (2k + 1) | - 1 ae 
9 T Ok(k + 1) 4k(k + 1)? 














50,7) =G2tb" —(2k + 1)r ul 
20,7) = “@" {er | 2k(k + 1) lle + ” 


= oy 2 —(2k + 1)n n 
+aeHl | 2k(k + 1) }Lare + 5] in}. (C.25) 


In particular, 2(0, ©) = 1. 
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BOOK REVIEWS 


Theoretical elasticity. By A. E. Green and W. Zerna. Oxford University Press, New 
York, 1954. xiii + 442 pp. $8.00. 


All serious students of the mathematical theory of elasticity will welcome this book. It does not make 
any claim to being a treatise on the subject. Rather, it represents both in choice of subject and presenta- 
tion the personal interests and point-of-view of the authors. Indeed many of the topics discussed are ones 
which are to a great extent the results of Professor Green’s own researches. Much of the material has not 
been previously presented in a unified notation or in book form. 

The subjects covered are drawn both from the finite and infinitesimal elasticity theories. In present- 
ing the basic concepts of these theories, finite elasticity theory is developed first, and the fundamental 
formulation of the infinitesimal theory is derived as a first approximation to the more general theory. 
Although this no doubt adds to the formal complexity of the presentation, it has the undoubted advan- 
tages of rigor and conceptual clarity. 

The finite elasticity theory is used to solve a number of simple problems which admit of exact general 
solution without further approximation or assumption—except that of incompressibility of the material. 
The theory of the superposition of infinitesimal deformations on finite deformations is then developed. 

The infinitesimal theory is developed for the cases of plane strain and generalised plane stress and is 
applied to the discussion of a wide range of problems for both isotropic and anisotropic materials. Complex 
variable methods are used consistently in a manner somewhat similar to that adopted by Muskhelishvili. 
In the final four chapters of the book the general bending theory of shells and membrane theory are 
briefly discussed, together with their particularisation to cylindrical shells and shells of revolution. 

Throughout the book tensor notation with convected coordinates is used. Although this makes for 
great elegance of presentation, it will no doubt limit to a large degree the extent to which the book will 


be used. 
R.S. Rrvirn 


Stochastic models for learning. By Robert R. Bush and Frederick Mosteller. John Wiley 
and Sons, Inc., New York, and Chapman & Hall, Ltd., London, 1955. xvi + 365 pp. 
$9.00. 


This book is primarily concerned with models for analysing learning experiments in which the learner 
is exposed repeatedly to the same stimulus. The essential ingredients of the model are relations between 
stimuli, responses and events. It is postulated that there is a set of probabilities associated with the 
various responses to the given stimulus. Following each response a certain event takes place such as 
for example a rewarding or an absence of a reward for the subject and this changes the response prob- 
abilities for future stimuli. The event will generally depend in some preassigned way upon the response, 
for example there may be a one to one correspondence between events and the responses. The events 
are represented by linear transformations to be applied to the response probabilities. 

Some very interesting mathematical problems arise in connection with the determination of the 
asymptotic behavior of the response probabilities after a large number of repetitions of the stimulus. 
Although the events themselves are represented by linear transformations on the probabilities, which 
event takes place depends upon the response to the previous stimulus; this in turn is a function of the 
response probabilities. The probabilities after any event are generally quadratic functions of those after 
the previous event 

Although the book contains many details that will be useful only to a psychologist who wants 
perhaps actually to apply the model to some data, it also contains numerous unsolved problems that 
make it very interesting reading for the non-specialist and it should greatly stimulate future study by 
both mathematicians and psychologists. 

G. F. NEWELL 
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SPHERICAL, CYLINDRICAL AND ONE-DIMENSIONAL GAS FLOWS* 
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JOSEPH B. KELLER 
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Abstract. The problem of spherical, cylindrical or planar flow of a polutropic gas 
has been formulated in Lagrangian variables, giving rise to a non-linear second order 
partial differential equation in two variables, A and ¢. A class of special solutions of this 
equation has been found by separation ‘of variables, and these solutions depend upon 
an arbitrary function which is related to the arbitrary entropy distribution in the gas. 
By specializing this function solutions corresponding to the expansion into a vacuum of 
an isentropic or non-isentropic gas cloud have been obtained as well as solutions cor- 
responding to the propagation of finite and of strong shocks in variable media. 

1. Introduction. Relatively few boundary value problems involving spherical or 
cylindrical flows of gases or compressible fluids have been solved exactly. In order to 
solve more problems of this type, particularly those involving variable entropy, we have 
investigated the flow differential equations, at first without regard to initial or boundary 
conditions. In this way a class of solutions of the differential equations, depending upon 
an arbitrary function, has been obtained. In general the solutions of this class represent 
flows with variable entropy, although some isentropic solutions are included. Then, the 
arbitrary function is adjusted to satisfy particular initial or boundary conditions. In 
this way the free expansion into a vacuum of a sphere, cylinder or slab of gas has been 
treated, as well as the propagation of finite and strong shocks in variable media. The 
latter treatment includes Primakoff’s point-blast solution as a special case. The expansion 
of a sphere of gas into a vacuum may be of interest to astrophysicists who have treated 
corresponding one-dimensional problems. 

The problem to be solved is formulated in Lagrangian variables, in Sec. 2 of this 
paper. In Sec. 3 a class of solutions is obtained by the method of separation of variables 
and in Sec. 4 the isentropic solutions are examined. In Sec. 5 the particle paths of the 
flows represented by the solutions are analyzed. In Sec. 6 the solutions are used to describe 
the expansion of a gas into a vacuum, in Sec. 7 to describe strong shock waves propagating 
in variable media and in Sec. 8 to describe finite shocks in variable media. Section 9 
contains conclusions and extensions of the results. 

2. Formulation. We consider the motion of an inviscid, non-heat-conducting fluid, 
obeying the polytropic equation of state. The one, two, and three-dimensional cases 
will be treated together. In the three-dimensional case y(h, ?) represents the radius at 
time ¢ of the particle with the Lagrangian coordinate h, which is defined by the equation’ 


wth,t) 
h= [ r"~"p(r, t) dr n= 1,2,3. (1) 
Jy(0.t) 
In the above equation p(r, ¢) is the density at time ¢ and radius r, and n is the dimension, 
which is 3 in the spherical case. In the two dimensional case y represents radial distance 
from an axis, and in one dimension y is a cartesian coordinate. 


*Received June 22, 1955. This work was sponsored by the Office of Naval Research under Contract 


No. Nonr-285(02). 
1R. Courant and K. O. Friedrichs, Supersonic flow and shock waves, Interscience Press, 1948, pp. 30-32. 
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It is fyrther assumed that all flow variables depend upon y and ¢ only and that flow 
occurs only in the y direction. This is the assumption of spherical, cylindrical or planar 


symmetry. 
From the definition of y, the particle velocity u is given by 


“u=¥Y,. (2) 

Similarly from (1) the density p and specific volume r are given by 
rT=p =y"'y. (3) 
Because of our assumptions concerning inviscidity and non-conduction, the entropy 
s of a particle is independent of time (at least between successive shocks). Thus we have 
s = s(h). (4) 
The function s(h) is given by initial data or by shock conditions, and is assumed to be 


known. 
The pressure p is given by the equation of state 


P = ple, 8) = g(r, 8). (5) 
For a polytropic gas or liquid 

g(t, 8) = go + A(s)r’. (6) 
The function A(s), the adiabatic exponent y, and the internal pressure g) are assumed 
to be known. 


In terms of the above defined quantities, the equation of motion is 
Yur = —y" [gy "Ysa + 9.81). (7) 
For a polytropic gas or liquid this becomes, using (6), 
yer = YA()(y"™"'ys) 7 "'(y""'ys)ny"”* — Arlyn) 7y"™. (8) 


Equation (8) is a second order partial differential equation for y(h, t). The coefficient A 
is assumed to be a known function of s, and by (4), of A. 

The problem we consider is that of finding solutions of (8). From a particular solu- 
tion, the flow variables can be found by using (2), (3), (4) and (5). 

3. Product solutions. Let us seek product solutions of (8) of the form 


y(h, ) = fri. (9) 

Inserting (9) into (8), and separating variables we obtain 
Per at (10) 
A IT Tr - Aw sy ST =>. (11) 


In these equations A is an arbitrary separation parameter and primes represent deriva- 
tives with respect to ¢ in (10) and with respect to h in (11). 
To solve (10) we multiply by 7’ and integrate, obtaining 


a 
UY = =n? 


2d logj +a (y = 1). (13) 


“M4 (#1) (12) 
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Here a is an integration constant. Thus, unless 7 is constant, which is possible only if 
\ = a = O or if j = 0, we have the solution 


fi lar ms a] dj=t (y¥#1) (14) 
ie Ln(1 — ¥) | 


| ertogj + ar"? aj = t (y = 1). (15) 


To solve (11) we differentiate out and obtain 
VAL P FM DSCSINL IY TT — POL SY TA! =. (16) 


We now consider the inverse function h = h(f) and denote h’(f) by g(f). Then (16) 
becomes 


yA[—Va ft tO — OP MUS TTT — LO ae TAD =a. (17) 
We shall first treat the case y ~ 1 by introducing z(f) and B(f) defined by 
z= q""', B(f) = A[h(f)). (18) 


With these definitions, (17) becomes (prime denoting differentiation with respect to f) 





2’ + | — — Iy - NF + (og B)’ | 


Y 
(19) 
Ay — 1) (m—1)(y—-1) 41 __ 
+ +B f —_ 0. 
The solution of (19) is (with G a constant) 
ms f 
p= somemnpe-nn| g — MY D f gprs ay |, (20) 
Thus we have g from (18) and (20), and finally since h’ = q, 
f _ f y/(y-1) 
h= / p'B"| G - Ay =D f —* ay | df. (21) 
fe 


Equation (21) gives f(h) implicitly, and thus provides a solution of (11) for y # 1. 
The solution may be written more simply by defining F(f) by 


FU) = E _ My =D [ jB? af”. (22) 
If \ ¥ 0 this can be solved for B(f) and yields 
Bf) = (-Af)"’")F. (23) 


The solution for the flow variables can now be computed from (2)-(5), (9), (21), (22). 
We have for y ~ 1 and \ ¥ 0, remembering that f = yj’ from (9), 


uly, t) = yj" (24) 
r(y, t) = —dyj"'/F'(yz’) (25) 
ry, t) = g +z ""F(yj). (26) 
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Equations (24)-(26) give the flow quantities. In fact, these expressions yield a solution 
of the Eulerian equations of motion for an arbitrary function F provided that j(t) is 
given by (14). This is the first main result of this paper. It is to be noted that in these 
solutions, u is proportional to y and that F must be so chosen that 7, given by (25), is 
positive. This implies that F must be monotonic in a region where y is of one sign. 

In the excluded case y ¥ 1, \ = 0 we have instead of (24)-(26) from (2)-(5), (9), and 
(21) 


u(y, t) = yj’j’ = yt, (27) 
y,) =fB’ yy CO’ = Co(yt”’), (28) 
WHy,O=go+7°°Q"''” =g+ lt”. (29) 


The new arbitrary function 6 and arbitrary constant / have been introduced and the 
expressions simplified by using (14), which gives 
jt) = a’”’t. (30) 
The origin of ¢ is chosen so that j7(0) = 0. 
The corresponding solutions when y = 1 are obtained from (17) by introducing 
B(f) as in (18), after which (17) becomes 


; ; ; 3 
q’'q '—(n— 1) f ao sa + " = . (31) 


The solution of (31) is 


f 


af) =f Bf) exp | — AfB°'(f) df. (32) 
Thus 
h(f) = | | eB) exp | — AfB'(f) ay | df. (33) 
Equation (33) yields the solution implicitly. Again we define 
F(f) = exp | — »fB™(f) df. (34) 
Then if \ + 0 we have 
B(f) = —dfF(F’)". (35) 


The solution for-the flow variables, ify = 1 and \ # 0, is again given by (24)-(26) with 
y = 1 and j(t) given by (15). Similarly for y = 1 and A = 0 the solution is given by 
(27)-(29) with y = 1, G = 1 and j(¢) given by (30). 

Equations (24)-(26) and (27)-(29) represent non-isentropic solutions of the flow 
equations depending upon an arbitrary function. In order to construct these solutions 
explicitly one need merely evaluate the integrals in (14) or (15). 

4. Isentropic case. The above solution can be specialized to the isentropic case by 
setting A = B = a constant. Then from (22) and (34) we obtain 


Ce, ( os \ Ss y¥/(y-1) 

F(f) = E ~ Ss LY f| y ~ 1, (36) 
Y 

F(f) = exp sue y =e, (37) 
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Thus the solutions given by (24)-(26) become, for y # 1 and A ¥ 0 


u(y, t) = yj’, (38) 
, n 1/ Y My oe 1) 2--2 eer 

rly, t) = 7'B ‘[o _ oy BUY YJ | ; (39) 
, *—ny Y r( _ 1) 9.9 hein leas 

p(y, t) = got) E ~ “a YJ ‘| . (40) 


In these solutions B and G are constants and 7 is given by (14). For \ = 0, and all y the 
solution (27)-(29) applies with b constant. 
For y = 1 and all A, (38) is unchanged, but the other equations become 


r(y, t) = 7"B exp a , (41) 
p(y, t) = go +7" exp Mt . (42) 


Here B is a constant and j is given by (15). 

As an example of these isentropic solutions, let 7(t) be given by (48) and assume 
G = 0 in (39) and (40). Then (38)-(40) yield the power solutions (if y # 1, n(y — 1)/2 
~—] 


9 





u(y, i) = - Se ap 
uy, t 7 D+ 2 ¥ P (43) 
B at? the 
—* E (oat 2)'] yey”, (44) 
a 9\2 }r/U-y 
ply, t) = go + p| (matt?) | ” & pate? (45) 


5. Particle paths. A particle designated by a fixed value of h, follows the path given 
by (9), ie. y = f(h)j(t) where j(t) is a solution of (12) if y ¥ 1 and of (13) ify = 1. 
Let us first consider the case in which the constant \ in (12) and (13) is zero. Then in 
both cases we have 


j(t) = 70) + at. (46) 


If the origin of ¢ is shifted so that j7(0) = 0 then this becomes 
y(t) = zat. (47) 


This solution has already been used to simplify (27)-(29), [see Eq. (30)]. In the motion 
described by (44) with the plus sign, all the particles start from the origin at ¢ = 0 and 
each moves outward with the constant speed a’’*f(h). With the minus sign all particles 
move toward the origin, each with the constant speed —a'’*f(h) and they all arrive 
there at ¢ = 0. 

Now let us examine the case y > 1, \ > 0. Then in order that (j’)’ in (12) be positive, 
we must have the constant a in (12) positive. From (12) we see that the solution is then 
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of the form shown in Fig. 1 when the origin of ¢ is chosen so that j’(0) = 0. This choice 
of origin requires that the constant a = —2Ajo°~”’/[n(1 — y)]. 

This solution is even in ¢. Each particle has the constant speed —a'f(h) at t = — © 
and slows down until it comes to momentary rest at t = 0 with y = f(h)j(0), after which 


t 
5 


i i i A. — 


0 1 2 3 “ t 


Fie. 1. The behavior of j(¢) for y > 1, \ > 0. The origin of ¢ is so chosen that j’(0) = 0 which requires 
@ = —2)j(0)""-” /[n(y — 1)]. Att ~+ &, j’ + ¥a!/*, The curve is obtained from Eq. (14) with y = 1.4, 
n= 3, = 7(0) = 1. 








it again moves outward ultimately approaching the speed +a’”f(h). The previous 
solution in which 7(0) = 0 is obviously a limiting form of this solution. If we only con- 
sider that half the solution ¢ > 0 we have the expansion of a gas initially at rest. 

The other possibility when y > 1 is that \ < 0. The nature of the solution now 
depends upon the sign of a as is shown in Fig. 2 which is obtained by analyzing (12). 
In the first case (a > 0) each particle has the constant speed —a’’’*f(h) att = — © and 
speeds up until it reaches the origin at ¢ = 0 with infinite speed. Thus all the gas col- 
lapses into the origin. The positive branch of the curve describes the reverse phenomenon 
—all the gas is released from the origin at £ = 0 with infinite speed and it spreads out 
and slows down, each particle approaching the constant speed a’”f(h) as t becomes 
infinite. The second case (a = 0) is similar to the first, but the speed approached as ¢ 
approaches plus or minus infinity is zero. The explicit solution in this case is 


J Dn 1/2 ny —- 1) +2 2/[n(y—1) +2) 
HD) ~ | = —_ | (" 2 )e 


ify ~ 1,a = 0, nq — 1) ¥ --1. (48) 


The third case (a < 0) is different from the preceding cases. The solution is periodic, 
all the gas being released from the origin at some instant, expanding, coming to rest, 
collapsing into the origin again, etc. 

When y = 1 and d < 0 the solution has, for all values of a, the behavior shown in 
Fig. 2c while for \ > 0 it has the behavior shown in Fig. 1, as can be seen from (13). 
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Figs. 2a and b. (For legend see p. 178.) 
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0 1.0 20 t 30 4.0 
Fic. 2c. The behavior of j(¢) for y > 1 and > 0. The origin of ¢ is so chosen that j(0) = 0. In the first 
case, a > 0,7’ approaches +a!/? ast ++ ~. In the second case a = 0, j’ approaches zero as t ~+ , and 


in the third case a < 0,7 is periodic in ¢. The curves are computed for y = 1.4,n = 3,A = —1,a = lin 
the first case and a = —1 in the third case, using Eq. (14). In the second case a = 0, the solution is 


given by (48). 


6. Application: Free expansion of gas into vacuum. Let us apply the isentropic 
solution (38-40) to the free expansion into a vacuum of a gas with y > 1. Then the pressure 
must be zero for some value of h, which represents the interface between gas and vacuum. 
If the interface is initially at y = ys , and if g. = O since the medium is a gas, we need 
merely set G = [A(y — 1)]y7(0)"?/2yB’’” in (39, 40). Then p = 0 at the interface which 
is given by 

et ee. ¢ 
Y = Yo (0) (49) 
The solution j(¢#) which must be taken is of the type shown in Fig. 1, if yo and 7(0) are 
not zero. Thus the gas is initially at rest and, according to (40), the initial pressure 
decreases from a maximum at the center y = 0 to zero at y = Yo as shown in Fig. 3. 
Then the gas accelerates outward, each particle approaching a constant speed while 
the speed at any instant is proportional to y. The pressure at any particle varies as the 
—ny power of its radius, thus decreasing as the radius increases. The preceding solution 
applies to the expansion into a vacuum of a sphere (n = 3), a cylinder (n = 2) or a slab 
(n = 1) of gas initially at rest, and with uniform entropy. 

In a similar manner, but by using solutions for j(t) of the types shown in Figs. 2a 
and 2b, solutions describing the expansion of an isentropic gas which fills all space 
except an expanding region of vacuum around the origin (i.e. a hole), can be constructed. 

For a gas with y = 1 which is also isentropic, it is possible to construct the solution 
for an expanding gas cloud which fills all space initially, but which has a Gaussian dis- 
tribution of density. To this end we choose \ > 0, in which case j(t) has the form shown 
in Fig. 1. Then wu, ¢, and p are given by (38, 41, 42), with g, = O in (42) if the medium 
is a gas. This solution describes a gas cloud initially at rest, which subsequently expands, 
the velocity of all particles ultimately approaching the same constant value. If this 
solution is also considered for t < 0 it describes a gas cloud which contracts, comes to 
rest and expands again. 

Another interesting motion of a gas with any value of y results if (30) is used for 
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7(t) with the plus sign, and if X = 0 in (39, 40). We then obtain, setting g. = 0 in (40) 


to describe a gas, 








u=yt', (50) 
r= a”’Br, (51) 
p=a”"B*¢"". (52) 

p 

1.0 

8 

6} 

4} 

ral 

Oo 5 1.0 5 20. y 25° 


Fic. 3. Initial pressure distribution in an expanding spherical, cylindrical or planar gas cloud, computed 
from Eq. (40) with g = 0, y = 1.4,G = B =j = 1, = 1. The same curve applies at later times if the 
vertical scale refers to pj”” and the horizontal scale to yj™. 


This solution represents the expansion (or contraction if ¢ is replaced by —?) of a gas of 
uniform density, pressure and entropy which fills space. 

Similar motions of gases in which the entropy is not initially constant can also be 
constructed by using equations (24-26) to describe the flow. In order that p equal zero 
at the interface between gas and vacuum, F must be chosen to vanish at some value of 
its argument, which value will then correspond to the interface. Aside from this condi- 
tion, F may be chosen arbitrarily. 

7. Application: Strong shocks in variable media. Let us suppose that a shock given 
by the equation y = R(t) moves into a variable medium of density p(y), pressure po(y) 
and velocity zero. The pressure p, velocity u and density p just behind the shock are 
related to the corresponding quantities in front of it by the shock conditions. These 


conditions are 











ep  (y¥+)p+0—-)Dp. y+1 . 
= Pr hy — Po YT 1 53 

es 2(p — Po) — | 2p | . (54) 
(2pol(y + Dp + Cy — Vpol}'” ~ La + Dood ’ 

lias | + bet (y - Dp. | a [ot Be | (55) 


The second expressions on the right apply to a strong shock, for which p > po . 
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We shall now assume that the flow behind the shock is a “product’’ solution given 
by (24)-(26), and that the shock is strong. We wish to determine the functions F(yj~"), 
R(t), po(y) and po(y) for which such a solution is possible. To this end, we insert (24)-(26) 
into (53)-(55) and obtain 











F(R7") _y+1 . 
—\Ri p(k) y—1 (56) 
Rj’ _ | 2 + af Fie | 57) 

j | ie 1) po(R) 

ae (y¥ + Ii go + pre) |” (5 

-* | Dpa(R) ; aad 
From (57) and (58) we have 

R (y¥+1) 7 = 

an ee ib (59) 


Solving, and introducing the integration constant Ry , we have for R(¢) the expression 
R(t) = Rf[x()]***””. (60) 


Inserting (60) into (57) and making use of (12), satisfied by 7(¢), we find 





2 — ¥) 


(+41) - 9/ _— 2r —I1,\ —2n 
F(z) a + 1 po(2 Y+1)/(¥ pp? (1 ”) _— (xRo 7 4 p| 
; 7 ha eal ae a tl eld — apr (61) 


Thus R is determined by (60), F is related to po by (61), and (57), (58) are satisfied, 
although the constants R, , D, and ) are still arbitrary. We now insert (60), (61) into 
(56) to determine F. We obtain, if g, = 0, 





a J 2,R?" a 1? 
F(a) = F x Fn a Dzx* | . (62) 
From (61) and (62) we then find an expression for p.(#), namely 
aa 7. 2r0Ro" (4n)/(y¥ +1) macrcanvevany 
ss E ro = 
P Jpte—O) y—-Sn ei Cy +t pp isin~8) year ei /ity otiy—s)) (63) 


In (62) and (63) F, is an arbitrary constant. The corresponding results with g, ~ 0 
are somewhat more complicated. 

We have thus obtained a solution with a strong shock moving into a variable medium 
at rest, with density given by (63). The constants F, , Ry, , \ and D are arbitrary in this 
equation, but only two essential combinations of these constants occur. The shock 
curve is given by (60), and j(t) by (14). The flow is given by (24)-(26) with F(z) given 
by (62). The flow might be produced by a piston following one of the particle paths 


Y(t) = Yoj(d). (64) 
The solution was deduced for y ~ +1 anddA = 0. 
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As an example of these solutions, let us suppose that p.(R) = constant. From (63) 
we find that this requires 





3n — 2 
D= 0, y= ca (65) 
The last condition can be fulfilled only for n = 3, in which case y = 7; (n = 1 leads to 
y = —1, which was excluded in the derivation, and n = 2 yields no value of y). In this 
case we have from (14), (62) 
_ 5(—v)'? ia =. il 3F Ro" 
a= | 3 t ’ F(z) =F (—n)'” (66) 


The solution computed by using (66) in (24)-(26) (with » < 0) is exactly the point 
blast solution of Taylor’s type first found by H. Primakoff, and applicable to high 
energy explosions in water. [See’, p. 424.] 

8. Application: Finite shocks in variable media. We may also apply the above 
method to determine finite shocks in variable media. Then we must satisfy the exact 
shock conditions (53)-(55) rather than the strong shock conditions. Inserting (24)-(26) 
into (53)-(55) yields with g, = 0 

k (x) a> Dy F + (y — 1po . saa (67) 
ART Po (y + Ippo + (y — GF 





Rj'j*" —_ 2: F ee Po (68) 


~~ (po/2)7 [ly + DIT"F + Gy — Vpol'?’ 





° 1 — 1/2 a 
R= Ga ler + DIF + Cy — Vpod'”. (69) 


Equations (67)-(69) are a set of two first order ordinary differential equations and 
one algebraic equation for the determination of the four functions F(z), R(t), po(R) and 
po(R). Equation (14) gives j(t). This system is evidently under-determined and will 
have infinitely many solutions. In the strong shock case, however, > did not occur and 
thus the system was determined. 

To solve the above system we could impose some relation between p) and p) and 
eliminate both these functions by means of that relation and (68). A pair of simultaneous 
first order equations for F(x) and R(t) would result. However, these equations can be 
treated separately since, from (68) and (69) we have by eliminating F, 


; Pa diana “— 7 
R = VES Ri (ore (Rj'j + rposs' | (70) 


This is an equation for R(t) if the ratio pop’ is given as a function of R. (Equation (70) 
also holds when g, ¥ 0.) After solving this, (67) can be solved for F, and then pp and po 
can be obtained. 

We will now investigate the shock curve R(t) given by (70), in the special case in 
which j(¢) is given by (48) with C = 0 and ypopo' = co is constant. The latter assumption 
means that the temperature ahead of the shock is constant. Making use of (48), (70) 
becomes 


- yt! 1 (y + 1)’ -n3 i 
B= eg ee + pt pary +e]. an 
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If we now introduce U(t) = Rt~ in (71) we obtain an equation for U in which the 
variables separate. One solution is the constant solution 


a de a ee. 5] Y 

= so ty~- +2 

In this case, which is physically possible only when (y + 1)/[n(y — 1) + 2] < 1, the 
shock curve is the straight line 

> y — . = = be a3 i sg 72 

R= Ut = +t (1 a= as t. (73) 


When U’ does not have the value given by (72), the equation for U yields 


+ F se ‘th +6] 5 | dU = log ; . (74) 


[2n(y — 1) + 4)? 


6) l 1#a—1 r) 1/(#a-a+l 
bt +a + (a 1) cos ¢ i »(sin 2) lit = (75) 


Here 6 = cot™'[aU/eo], a (y + 1)/[2n(y — 1) + 4] and bisa constant. Thus we have 


t 


R= Ut =. bt cot 8. (76) 
ab 
Since bi is given in terms of @ by (75), R(t) is given parametrically in terms of 6. A graph 
of [(ab),c.]R versus bt is given in Fig. 4 for n = 3 and y = 1.4, which represents an 
expanding spherical shock in air. When the minus sign is chosen in (71)-(75) the same 
curves are obtained as with the plus sign, provided ¢ is replaced by minus ¢, thus repre- 
senting a converging shock. Of course, only one of any such pair of solutions can occur 
physically, since only one obeys the entropy inequality at the shock. 
For small values of ¢ we have 
R ~ Kt" (77) 


where K is a constant. For y = 7 and n = 3 we find a = 1/5 and thus R ~ ¢”’, which 
is the behavior of Primakoff’s point blast solution for smal] ¢. For large t, R behaves 
linearly in ¢, as in (73). 


[To determine p, and F, we have from (69) and the condition ypop. = Co , 
y Sig We a 1 
7 ap R cael '« Mieeaes l Do “yp °2 ony > dence ny ao 
F p I = ae YPo —R j an BB DoJ a (78) 
(y+ 1) °"’ y + Ie + I 


Thus, if dot denotes time derivative and prime denotes derivative of a function with 
respect to its argument, we have 


-_ " ; 
i rg 3 (ppk“j"* + 2p.R'R'j"” + nypR°S ‘’) 
) Ne? 
I 


F(R j"' — Rj-?j’). (79) 
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Fic. 4. The radius R of an expanding spherical shock in air as a function of time t computed from the 


parametric equations (75, 76) with n = 3, y = 1.4. 
Using (78) and (79) in (67) yields the following equation for po 


9 
__ “7 — 1°37 ‘ *D** mY *2my-l-y 
E +E De Poh I” + 2pok Ri" + nypokt PF’) 


\ 2 
)Co 


5 Pome ] atta? ey my—liy m= i 2 ° 1 “2 9 - 
sas —— (poRj"” + nypoj”'j ) tai ype (Rj — Rij’)T' 





9...-2p*2 coe 1)2yR” = 1)? |" > 
= 2y05°R*| y +14 G + De - a0" . (80) 
Upon introducing p, = piR’, the above equation can be integrated yielding p,[R()] 
as a function of ¢. To obtain p,.(R), we must insert ¢ = ¢(R) from (75) and (76) into this 
result. Then when p,[R(é)] is known, F can be found from (78). These calculations will 
not be carried out here. 

9. Conclusion and extensions. The problem of spherical, cylindrical or planar flow 
of a polytropic gas has been formulated in Lagrangian variables, giving rise to a non- 
linear second order partial differential equation in two variables, h and ¢. A class of 
special solutions of this equation has been found by separation of variables, and these 
solutions depend upon an arbitrary function which is related to the arbitrary entropy 
distribution in the gas. By specializing this function solutions corresponding to the 
expansion into a vacuum of an isentropic or non-isentropic gas cloud have been obtained 
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as well as solutions corresponding to the propagation of finite and of strong shocks in 


variable media. 
Other types of gas motion can be described by solutions of other forms. For example, 


solutions of the form y = f(at + Bh), where a and 8 are constants, can easily be found, 
and they describe steady motions. Solutions of the progressing wave type, similar to 
those considered by G. Guderley, by G. I. Taylor, by J. von Neumann and J. Calkin, 
and by R. Courant and K. O. Friedrichs (loc. cit., Chapter VI-C), can be obtained if 


y has the form 
y = tf[K(hjt]. (81) 
Here a is a constant and K(h) is a function of h to be determined. By inserting (81) into 
(8) we find that K(h) must have the form 
K(h) = (Bh+ 0)”°°", (82) 


where B, C and « are constants. 
The resulting equation for f can be reduced to first order, following von Neumann 


and Calkin, by introducing the new variables s, F(s) and q(F) defined by 


s = log [K(A)t], 
F(s) = e~**f(e'), (83) 
q(F) = (4 4 s)FO). 


In (83), the constant 6 = [an(1 — y) — 2a + 2][2 + n(y — 1)]~’. The equation (8) 
for f now becomes the following equation for q(F) 

[L — yDFOP Pg Yq — Fg’ 

= FPO” Dive + y(n — 1)QF"* + ey +a4+1—-—€] — (a — LaF — 2aqg. (84) 
In (84) D is a constant. Special cases of (84) have been studied by von Neumann and 
Calkin. 

Note added in proof. The solution (24)-(26) for y ¥ 1 and \ # 0, when specialized 
to the case g, = 0, coincides with the result of L. I. Sedov, ‘On the integration of the 
equations of one-dimensional gas motion’, Doklady, AN USSR 90, 735 (1953). Some 
related results are given by M. L. Lidov, “Exact solution of the equations of one-dimen- 
sional unsteady gas motion taking into account Newtonian gravitational forces”, Doklady, 
AN USSR 97, 409-410 (1954). 
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EXPRESSION OF WAVE FUNCTIONS OVER A HALF SPACE IN 
TERMS OF THEIR BOUNDARY VALUES* 


BY 
H. PORITSKY 
General Electric Company, Schenectady, New York 


1. Introduction. We consider the expression of a solution u of the wave equation 


ee. 
Vu — 35 = 0 (1.1) 


over the half space z > 0, in terms of its boundary values over the boundary z = 0: 


1 1 | du 
u(P, t) = ule, y,2,) = 5 il fu] do + = | fe aw, (1.2) 
and the alternative expression of u over z > 0 in terms of the boundary values of (du/dz) 
given by: 
— —1 [| | dz’ dy’, 
0 On [2 R ad) 


Here the integration is carried out over the plane z = 0, dw denoting the element of 
solid angle subtended at P by the plane element dz’dy’ at P’, while 


[u], [du/dt], [du/dz] (1.4) 


denote the retarded values of these functions at P’, that is their values at P’ not at the 
time t, but at the time t’ = t — R/c, where R is the distance PP’. The “retarded time”’ 


t’=t—R/c (1.5) 


R 








fK 





P 
(dx dy , 

! 

! 


wT 


Fic. 1. 





is the time ¢’ such that a spherical wavelet, starting from P’ at ¢’ and with a radius 
expanding with a velocity c, would reach the point P at the time ¢, at which the left 
member of Eq. (1.2) is evaluated. 
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In the version of this paper first submitted to this journal, the authors derived 
several proofs of these equations, one based on Fourier integral time-resolution, another 
on a slight modification of the method used by Kirchhoff in expressing solutions of the 


non-homogeneous wave equation 
Vu- 2a * —Arp(z, y, z, t) (1.6) 
“source” distribution p, as a sum of the retarded potentials of the latter, 


in terms of the 
as follows (see for instance [1]): 


“u= le] dv. (1.7) 

After this paper was revised in accordance with the recommendations of the referee, 
it was pointed out to the authors that Eq. (1.2), (1.3) appear in reference [2], where a 
hint regarding their proof is also given. It, therefore, seemed pertinent to omit proofs 
of (1.2), (1.3), and confine ourselves to several applications of these equations and their 
n-dimensional counterparts. 

The boundary value expressions under discussion, it should be pointed out, differ 
from the more commonly considered solution of the Cauchy problem for (1.1), in which 
u is expressed for t > 0 in terms of the values of u, du/dt at the time ¢ = 0. 

It is to be emphasized that the expressions S, of the wave functions in terms of their 
boundary values considered here are unique only if, as is customary in physical appli- 
cations, the retarded potentials are used. A similar, but different expression S, , utilizing 
the advanced potentials also exists. Thus a solution of (1.6) exists of the form (1.7) but with 
[o]/R replaced by the advanced potential p(t + R/c)/R. A linear combination of these solu- 
tions of the form AS, + (1 — A)S, for any constant d also furnishes a possible solution. 

In the special case where u is independent of time, Eq. (1.1) reduces to the Laplace 


equation, and Eqs. (1.2), (1.3) simplify to 





u(P, t) = = / u dw, (1.8) 
uP, ) = -~ [Sew (1.9) 


Both of these equations are familiar from potential theory. 
2. Examples. A. Azially symmetric product solutions of the wave equation. As a 
first example we consider the axially symmetric wave function 
u = J,(dr) exp [—(7y’ — k*)'’z] exp (iat), 3 r=a2rty’, (2.1) 
for y > k. Applying Eq. (1.3) for a point P on the positive z-axis, there results, upon 
introduction of polar coordinates in the plane of integration and cancellation of the 


factor exp (zat), 








plo eal. [ seehaor de (2.2) 
(y? a kK)? ‘ R 


Since (see Fig. 2) 
R=r'+2, RdR=rdr, (2.3) 
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one obtains 











ep LO tal = [exp (ith JutuiR? — 2)")dR. (24) 
R 
z 
a 
Fia. 2. 


In a similar way one obtains for the same wave function (2.1), by applying either 
Eq. (1.2), the following integral 


exp [—(7? — k)"”2] [ Jo(yr) exp (—ikR)(1 + ikR)erR™ dr 


(2.5) 


Z i exp (—ikR)(1 + tkR)Jo[y(R* — 2°)'”*] dR/R’. 


The integral (2.2) is essentially equivalent to (see Sommerfeld in [3], Watson [4]) 
exp [ik(y* + 2*)'”] = r wn f—(.2 — B2\/251/,2 2) —1/2 
LP aye — = | ddJolda) exp [-O" — H)72]0" — KYA. (2.6) 


Indeed, by carrying out the following substitutions in (2.6): 





ik — —z, A> Yr, r—y¥, z— tk, (2.7) 


one transforms (2.6) into (2.5). 

Further integrals can be obtained from the wave function (2.1) by choosing P in 
Eqs. (1.2), (1.3) not on the axis of symmetry. These integrals remain two-dimensional 
since the integrand is not axially symmetric about a line through P normal to the plane 
z = 0. However, by the application of the Neuman expansion the integral can be reduced 
to the axially symmetric case. 

B. General axially symmetric wave functions. As our second example we consider 
any wave function u(r, z, t) in z > 0 which is axially symmetric about the z-axis, that is, 


a solution of 


i) 
= 
Q 
= 


i 
r 


2 1 Fs 
3 -3 32 « @. (2.8) 


a 


= 1 
| 


eu 
ort 


Q 
& 
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Eq. (1.3) yields, upon replacing r by (R* — z*)’”, the following expression of u for a 
point on the positive z-axis in terms of u, = du/dz along z = 0: 


u(0, 2, t) = -f u,[(R? — 2*)”?, 0, t — R/c] dR. (2.9) 
Eq. (1.2) yields similarly 


w0,s,0 = z | ul(R? — 2)”, 0, t — R/cJR™ dR 
. (2.10) 


+ (z/c) [ u[(R? — 2*)'”, 0, t — R/c] dR. 


C. Plane waves. As a third example consider an arbitrary plane wave solution of the 
form 


u = f(t — 2/c) (2.11) 
where f satisfies the conditions 
f(s) = 0 for s <0, (2.12) 
f'® =0. 


Since the field is plane, it can be considered to be symmetric about the z-axis. Application 


of Eq. (1.3) for 
t—2z/ce>0, z>0 (2.13) 


yields, upon introducing polar coordinates as in examples A and B, 


| 


flit ~ s/) = —(1/6 [ ” $(t — R/orR™ dr | 
0 (2.14 


—(1e) [st — R/o aR, 


where 
tr as [(ct)? = oT. 


That this is an identity follows immediately upon carrying out the integration. 
Similar application of Eq. (1.2) to the plane wave function (2.11) satisfying (2.12) 
yields for z, ¢ satisfying (2.13) 


f(t — z/ce) =2 / . [f(t — R/)R™* + f(t — R/de'R™)r dr, 
, (2.15) 


- f f(t — R/OR™ dR + 207 [ f(t — R/R™ aR, 


where 7; is as in (2.14). That this again is an identity follows by integrating the last 


integral by parts and utilizing (2.12). 
It will be noted that the wave function (2.11) does not include a simple standing 


plane wave such as 
u = sin at sin kz. (2.16) 


Indeed for this example the right-hand member of Eq. (1.2) vanishes and Eq. (1.2) 
obviously fails to hold. The function (2.16), if resolved into travelling waves, includes 
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a plane wave travelling in direction of negative z. The reader will discover in supplying 
the proof of theorem (1.2) that the integrals over bl a hemisphere of large radius R; , 
which must converge to zero as Rz — @ to yield (1.2), fail to do so for the plane wave 
travelling in the direction of negative z. 

D. Spherical wave. As our next example we consider a spherical wave originating 


from a point P, ,z = — a,a > 0, on the negative z-axis (see Fig. 3): 
— f(t — fet ; (2.17) 
2 


where R, is the distance from P, . Again we assume that f satisfies the conditions (2.12). 


Pp 
we 











i 
P, 
Fie. 3. 

Eq. (1.2) now yields for a point P at z = b on the positive z-axis for 
t— (a+ b)/e>0, (2.18) 
upon recalling Eq. (2.10), the following relation 


fit — a+ Die] _ yf” Lh, ER] aR 
(a + b) 6 R,R* 














(2.19) 
" ae ft — (Re + R)/c] dR 
c Js RR 
The upper limit R’ of the integrations is given by 
R’ + R, = ct. (2.20) 
Upon utilizing the relations 
R=+r, R=a+r (2.21) 
one obtains from (2.19) 
= R,’ 
Lt— a+ De} sf” ste — a + R/S — (o? + BDI” aR, 
iii (2.22) 


R,’ 


—(o/) | f'lé- a+ R)/c)[R2 — (a’ + b’)] dR, , 
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where 
Ri = (a — b’)/2ct + ct/2. (2.23) 


That (2.22) is an identity follows by integrating the second integral by parts and utilizing 
the initial conditions (2.12). 

3. Extension to n-dimensions. We now consider the extension of the above results 
to n-dimensional (Euclidean) spaces E, and to solutions of the ‘“‘wave equation” 





au eu ’u 1 du 
asataat::: 3-357 =0 3.1 
Oz; * Ox; + +r i nn ae) (3.1) 
The aim is to express u over the region z, > 0 in terms of its boundary values over 
z, = 0. 
The situation is quite different depending upon whether n is odd or even. For odd 
n an extension of Eq. (1.2) 7s possible and u(P, ¢) can be expressed in terms of retarded 


values of u, du/dt, --- , 9°” u/at®*”” along the flat z, = 0 as follows: 
9 e (n-1)/2 C _ on"u 
uP, ) == | want | SS | du, « 3.2 
Q,, : p> c™ ot” Ww ( ) 
Here C,,, ; m = 0, 1, --- , are proper constants, the integration is carried out along 


xz, = 0, dw, denoting the element of “‘solid angle’ subtended by dz, --- dz,_, at P’ 
in z, = 0 at the point P: 


di, = za dxi +++ dat, , (3.3) 
while [u], [du/dt] , --- , as for n = 3, are the values of u at P’ at the time ¢ — R/c, 
where R = PP’; finally 
a [ = n{T(1/2)'/T[(n/2) + 1] (3.4) 
is the complete solid angle, that is the (n — 1)-dimensional content of the “‘surface’”’ 
of the n-dimensional unit sphere: 
et+at ee te = 1. (3.5) 


Similarly, for odd n, a generalization of Eq. (1.3) is possible, and is given by 


9 (n—3)/2 Jom am+1 
u(P, t) = os eee / p> ee ES dx} a he dx} -, ’ (3.6) 
where again D,,, are proper constants. 

For even n the situation is quite different, and the expression of u in terms of its 
boundary values is more complicated. This is but one instance of the many differences 
between solutions of the wave equation in even and odd number of dimensions. In the 
present case, one may ascribe this difference to the fact that the Bessel functions which 
are involved in the proof of Eqs. (3.2) or (3.6) turn out to be of integral order for even 
n, but of an order that differs from an integer by a half for odd n. The latter Bessel 
functions, as is well known, can be expressed in terms of exponentials and a product of 
a finite number of negative fractional powers of the argument; no such simple expression, 
however, is possible for the Bessel functions of integral order. 


1956] WAVE FUNCTIONS OVER A HALF SPACE 191 


The expression of u in terms of its boundary values for even n is considered briefly 
toward the end of this section, and in more detail in the following section. 

After these preliminaries we indicate briefly the generalizations of (1.2), (1.3) for any 
number of dimensions n by a method similar to the Fourier integral method mentioned 
in Sec. 1. We resolve u, the solution of (3.1), into a Fourier integral in time: 


u= [exp (at)U(er , 22, +++, 2) day, (3.7) 
where (it is assumed that) the integrand satisfies Eq. (3.1), and hence U is a solution 
of the n-dimensional equation 

V7U + kU =0, = a/c. (3.8) 
We apply Green’s theorem to U and the Green’s function G for the half space x, > 0 
and the condition U = 0 onz, = 0: 


GP, P’) = falk, R) — falk, Ri), (3.9) 


where F is the distance from P, R, the distance from its mirror image in z, = 0, while 
f,(k, R) is a proper, spherically symmetric solution of (3.8), which behaves like 


Const. exp (ikR)/R°"”” (3.10) 
at infinity and like 
2-n =_ 
I" /n — 2) for n> 2, 3.11) 
—lnR for n = 2, 


near R = 0. There results 


2 


UP) = = foe f xaU ery tay +++ te, OL AfaCk, R)/R AR) det +++ dat , (8.12) 
and, upon utilization of Eq. (2.3), 
U(P) = = f _ il Ulay , 22, +** , tar, OR" Of.(k, R)/@R deo, . (3.13) 


In a similar way, by replacing G by the sum of the two right-hand terms in Eq. 
(3.9), one obtains 


U(P) = 2/ i | U.Ats » ay °** Zana, Ofek, R) dt +++ dats. (3.14) 


To obtain f, we note that for spherically symmetric solutions Eq. (3.8) reduces to 


aU ,n—19U 
oR’ R OR 








+ kU =0. (3.15) 


Solutions of this equation are given by 
U = RB, 2-1(kR) (3.16) 
where B, is a Bessel function of order ». We shall choose 
f. = Const R'*?HY2_,(kR), (3.17) 
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since the asymptotic form of this solution yields a behavior at infinity of the form (3.10) 
which is more in accord with a “retarded potential’’ solution corresponding to an ex- 
panding spherical wave, than would be the case with any other Bessel function. 

The difference between even and odd n appears when we examine the expansions 


for H,,._, (see [3], p. 198, Eq. (6)): 





2\ : yA 1/2 | - 
~~ (3.18) 
(v, m) = (4r° — 1°)(4v? — 3) +++ [40° — (2m — 17] 
| 2°"m! 
=n/2—1 


For general v these expansions are asymptotic, and divergent. However, if v 
where 7 is an odd integer, then (v, m) = 0 for m 2 (n — 1)/2, the series in (3.18) reduces 
to (n — 1)/2 terms, and yields an exact representation for H<”’. 

Hence for odd n, one obtains from Eqs. (3.16), (3.17), upon adjusting the multi- 
plicative constant in accordance with Eq. (3.11), 


(n-3)/2 


R* exp (—7kR) >. ((n/2) = ‘ m)(2ikR)*~-?”?-™ 
m=0 


fk, R) = G—pyim/2) — 1, (n — 3/2] 





(3.19) 
(n—3)/2 


= exp (ikR) R?" >> Din(tkR)”. 
m=0 


Substituting this expansion in Eq. (3.13) and recalling the factor exp (tat) one obtains 


U(P) exp (iat) 
9p e (n—3)/2 (3.20) 
= -3 | tae | Uist, Me. ** - tla, 8 p> exp [i(at — kR)]D,,(ikR)”. 


Now replace the left-hand member by u(P, ¢); note that multiplication by the factor 
1k is equivalent to application of the operator (1/c)(@/dt), and that ¢ occurs only in the 
form exp [ia(t — R/c)]. Hence for functions u of the form exp (tat) U(x, --- z,) Eq. 


(3.20) may be recast in the form (3.6). 
For general wave functions, Eq. (3.6) now follows by means of Fourier integral 


superposition, as in Eq. (3.7). 
A similar proof of Eq. (3.2) follows from Eq. (3.12). By differentiation of (3.19) 


one obtains 


a) (n-—1)/2 ; 
9fn — exp (ik RVR" Y Cun(ikR)", (3.21) 
OR m=0 
where 
be 5 = 2 — 0 
= ( n) D, P (3.22) 
Cam = kD, m-1 + (m — n)Dya.n ; m> 0. 


The case n = 3 has been considered in Sec. 1. For this case the series in Eq. (3.18) 
reduces to its first term, unity, and Eq. (3.19) yields 
fs = exp (—7kR)/R. (3.23) 
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For n = 5, 7 one obtains 
folk, R) = SPE ty + ane, (3.24) 
= +7. ° 2 
tik, « Se aoe [1 + ue rm ae | (3.25) 


It is of interest to note that for odd n the functions (3.19), multiplied by the factor 
exp (zat), can be written in the form 





b (n—3)/2 R™D 3” ( R) 
on 2-n ae a 9 
u(R,t) =R p> ” oP F\t (3.26) 
where 
F(u) = exp (tau). (3.27) 


Hence, by superposition over a, it follows that (3.26) is a solution of the wave equation 
(3.1) for an arbitrary function F. This is the general expanding spherical wave, and is 
analogous to F(t — R/c)/R for n = 3, to which, indeed, it reduces for n = 3. 

Returning to the case of even n, now one can no longer utilize the asymptotic Eq. 
(3.18) to advantage for the integral order Bessel functions. Equations (3.17), (3.11) 
yield (for both even and odd n) 


miH,(aR/c) for n = 2, 
(3.28) 





f = a (n/2—1) eee 
Jn m (2 Re? Hn-1(aR/c) for n> 2. 


I (n — a(2 os 2). 


Again, Eqs. (3.13), (3.14) can be used, combined with Fourier a-integration as in Eq. 
(3.7), to yield the desired expression of u in terms of its boundary values, but the result 
is rather complex. More attractive forms are considered in the following section. 

4. Alternative methods in n dimensions. With the odd-dimensional case disposed 
of, one may treat the case of even n by immersing the space FE, of (xz, , %2 , *** , %n) in 
an odd-dimension space E,,, obtained by adding one coordinate z,,, to Z, , and regarding 
the solution u of Eq. (3.1) as a solution in E,,, of the corresponding (n + 1)-dimensional 
wave equation. This is sometimes known as the “principle of descent” (see [5]). 

As a first application of this “principle” we shall obtain for even n a general ex- 
pression for u,(R, t), and expanding spherical wave. This will be done by immersing 
E,, in E,,, , lining the z,,,-axis with a uniform distribution per unit z,,, of (n + 1)- 
dimensional point “sources”, and adding the resulting wave functions given by Eq. 
(3.26) with n replaced by (n + 1), for the elements of the wave components. Thus, for 
n = 2, we immerse the (z, y)-plane in an (z, y, z)-space, and obtain a wave of cylindrical 
symmetry by adding spherical wave elements originating from sources distributed 
uniformly over the z-axis. 

For general even n the resulting wave function is independent of z,,, and can be 
evaluated in the flat z,,, = 0. There results 


+2 n/2-1 


u,(R, t) = 2/ (RO Dyas el F(t — RM /0) dinar , (4.1) 


—-2 m=0 


where 
R” _— [R? + eo a (4.2) 
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Changing to s = R” asa variable of integration, one obtains 
t = ta] 
i) 2+m—n 
8 D2 ) 
u(R, i) = 2 5 ___a+im p™)(¢ — g/e) ds: (4.3) 
2 oe — BR) ; 
while introduction of 
: R sinh v, s = R” = RK cosh», (4.4) 
yields 
u(R, t) = >> Dasi.m(R cosh v)?*"-"c"F™ [t — (R/c) cosh v] dv. (4.5) 


m 


2, Fig. 4 shows in perspective the (x, y)-plane, its “immersion” in the 


For n = 
and the various distance R, R’’, dx; = dz involved in Eqs. (4.1), (4.2). 


3-space E; , 








te 
dz 2 
/\ | 
/ ~ 
x/ 
/ 
Fic. 4 
If F satisfies the conditions 
$s) = or s< 
F(s) 0 for s< 0, (4.6) 
f{® = F’(0 =--- F*’'"O), 


the upper limits in the integrations in Eqs. (4.3), (4.5) are replaced by ct, cosh™’ (ct/R). 
For such a case, for n = 2, Fig. 5 shows for a fixed ¢, both in perspective and in the 
plane z = 0, the spherical wave fronts corresponding to the argument 0 for F. 

For n = 2, 4, (with (4.6) holding), Eqs. (4.1)-(4.5) yield 


ct 
u(R, t) = 2 [ F(t — s/c)(s’ — R’*)~*’*ds 
a (4.7) 


osh t/R 
= 2 [ F[t — (R/c) cosh v] dv; 


J 


d 
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9 r*'! F(t — s/e "(t — s/c f 
u(R,t) == ll E 2 xe ~ tS — 2/0 73 _— 
3 JR § 8c (s Rk ) (4.8) 


= EE — (R/c) cosh») _ RF'(t — (R/c) cosh | d 
a A Sioa. so v. 
cosh” v cosh 





bo 





It is possible by integration by parts to eliminate the derivatives F’, F’’, from Eqs. 


(4.3), (4.5) leaving only F in the integrand. 











It will be noticed from Fig. 5 and Eqs. (4.1)-(4.8) that if F satisfies the conditions 
(3.2), the wave function u,(R, t) does not vanish outside the spherical shell 


0<R/e-t< i. (4.9) 


In this respect the case of even n differs from the odd one. 

Of special interest, both for odd and even n, is the limiting wave function V,(R, ¢) 
approached by w,(R, t) by letting F above approach a proper constant multiple of the 
Dirac 6-function: 

F(t) > C&(t). (4.10) 
The function V,(R, #) may be regarded as the retarded potential solution due to pulse 
point source at the origin, at time ¢ = 0, that is, it is the “‘retarded” solution of 


‘. 9° ' ; 
Vu-— uy i = —(,5(2, ,%2,°°* », 2,28); (4.11) 
c ot 
where 
5(2, es t) = 6(2;) 6(22) wth as 6(t). (4.12) 


The functions V,, for n even or odd, are useful in deriving the following n-dimensional 
analogue of Kirchhoff’s theorem. The solution of wave equation 
au au 1 du 
atte + es = -—D,p(t, , «°° » te, fe 4.13 
Ox; + Ox; Cc at nP( 1» »>“n ) ) ( ) 
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with a “source” P, is given by 
u(P,t) = [at | aP’ p(P’, t)V(R, t — ¢) dP’ (4.14) 


where 
R = PP’. (4.15) 


For odd 7, recalling Eqs. (3.26), (4.10) and replacing F(t) by 6(¢), one may convert the 
right-hand member of Eq. (4.14) to a form resembling Eq. (1.7); for even n, however, 
the existence of the ‘‘tail’” in the wave element originating from each source element 
prevents a true retarded potential solution from ever acquiring a form similar to Eq. 
(1.7). 
Finally, the function V,(R, ¢) can be used to obtain analogues of Eqs. (1.2), (1.3) 
for n-dimensions. This will be illustrated for n = 2. Equations (4.7) now yield 


2 2-1/2 
f (ct) — R’] for Rk < cet, (4.16) 
0 for R < ct. 


V.(R, t) = 


By utilizing this function, one obtains the following expression in y > 0 of the solution 
u of 
42 a 2 
Ie a (4.17) 
Ox Oy 
in terms of its values along y = 0, along the lines of Eqs. (1.2), (1.3): 


u(x, y) = : [ av [ dx’{u(x’, 0, t’)y V(R, t — t’)/R OR}, (4.18) 


5 T 


2 [ av [ dx’'[u,(x’, 0, t’)Vx(R, t — t’)], (4.19) 


where 
R?=(¢-27'//+y’. 


An alternative way of arriving at the expression of the two-dimensional wave function 
u(z, y) in y > O in terms of its values on y = 0, is again by applying the ‘‘method of 
descent”’, by considering u as a wave function in three dimensions, and applying Eqs. 
(1.1), (1.3) to the half space y > 0. First one obtains, since u is independent of z, 


u(x, y, t) = (1/27) [| [u(x’, 0, t — R’’/c) y/(R’’) 


+ u(x’, 0, t — R’’/c) y/cR’’) dx’ dz’, (4.20) 


u(z, y, ) = -3- [/ [u,(x’, 0, ¢ — R’’/c)/R’’) dx’ dz’ 
where 
(R"P =(¢4-—27yPt+yt()’. (4.21) 
By a proper change of variables (and integration by parts) it is possible to convert the 
2’-integration of Eqs. (4.20), (4.21) into a ¢’-integration as in Eqs. (4.18), (4.19), and 
derive the latter equations from the former. 
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—NOTES — 


CONCERNING THE EIGENVALUES OF A 
DIFFERENTIAL EQUATION IN CONVECTIVE HEAT TRANSFER* 


By HENRY E. FETTIS 


(Aeronautical Research Laboratory, Wright-Patterson Air Force Base, Ohio) 
In a recent paper by Batchelor [1] there occurs the differential equation 
d’f/dy’ + d»yG y)f 0. (1) 


It is required to determine the values of \ for which (1) possesses a non-trivial solution, 
subject to the conditions 


f(0) f(4) 0. (2) 


An approximate value of A was obtained in the above cited paper by expanding “f” 
in & power series 1n tniortun itely, ut h an expansion is too slowly convergent to rive 
a good estimate of A with only a limited number of terms. However, if the function f 
is expanded as a power series in the eigenvalue A, a much better value can be obtained 
with an equivalent number of terms. The means of obtaining such a series is quite simple 
and in fact a well known method, although its application to eigenvalue problems does 
not appear to have received the attention which would seem warranted in view of its 


utility as a means of obtaining numerical results. 
A more convenient form of the equation for this method is obtained by setting 


y (¢ + 1)/4. Denoting derivatives with respect to “t’’ by primes, the new equation is 
{"+a(l O)f = 0, (3) 
where a \/(16)°, and the boundary conditions are f(+ 1) = 0. It may be noted that 


the general solution of Eq. (3) may be given in terms of the confluent hypergeometric 


function [2]. However, because of the limited tabulation of this function, this form of 


the solution is of little value in the present eigenvalue problem. Therefore, proceeding 

according to the previously mentioned plan of obtaining a formal expansion as a power 
series in a’, we assume that f(t) may be written in the form 

P . , le ’ in \ 

f fit) taf) taflt) +--+ +. (4) 


If this series is inserted in the differential equation (3), and the coefficients of the various 
powers of a’ set equal to zero, there results the following series of equations: 


T° 0. Bees (?? — LDYo% rete) (? — lfi, etc. (5) 


Solving this system in succession gives the desired coefficients in (8) as functions of “‘?’’. 
The boundary conditions require that f(1) 0, which results in the characteristic equa- 


tion in the form 


fo(1) + a’ f (1) + a f.(1) + ++: + = 0. (6) 


*Received Feb. 21, 1955; Revised manuscript received March 30, 1955. 
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If the series is terminated with the nth term, then (6) is merely a polynomial equation 
in a whose roots are approximates to the eigenvalues. 
The solutions of (5) for the case where f is even in ¢ are given below 


Sold l, 
f(t) = (—.5? + .8333333¢') x 10°", 
f(t) = (4.1666667* — 1.94444440° + .14880950") x 10°’, 
f(t) = (—1.3888889¢° + 1.09126970° — .23258372'° + .01127351'*) x 10™, 
f(t) = (.2480158¢° — .2755732t'° + .1002917t'? — .01339877'* + .0004697t"*) x 10~, 
f(t) = (—.275573t"° 4+ 04123671’? — .0206519¢"* + .00473711"° 
— 000415322" + .00001242"") « 107°, 

fo(t) = (.0020877t'? — .0037799t'* + .0025787'° — .00082977"* 

4+ 00013662" — .00001017°? + .0000000277*) « 107°. 


Hence 


fo(1) 1, 
fi) 1.166667 X 10°', f,(1) = 2.3710318 X 1077, f,(1) = —.5189294 x 10°, 
f.(1) = 0598053 X 10, —f,(1) = —.0042472 x 107°, —sfa(1) = .0002044 x 10°. 
If only terms through a® are retained the equation 

8° — 4.16666678* + 2.37103188 — .5189294 = 0 (7) 
is obtained, where for convenience we have set 10/a° = 8. The largest root of (7) is 


8, = 3.528, 
and the corresponding values of a’ and ) are 
a; = 2.834, A, = 725.6. 
Listed below are the corresponding results obtained when 5, 6, and 7 terms are carried 


in the series: 


Five terms: 


8° — 4.166666768° + 2.37103188* — .51892948 + .0528053 = 0, 
8B, = 3.53638, a; = 2.82775, A, = 723.9. (8) 
six terms: 
8° — 4.16666678* + 2.37103186" — .51892948" + .05980538 — .0042472 = 0, 
B, = 3.536265, a; = 2.827762, x, = 723.91. (9) 


Seven terms: 
8° — 4.16666678° + 2.37103186* — .5189294," 
+ .05980538 — .0042472 + .0002046 = 0, 
8, = 3.5363645, a; = 2.827763, A, = 723.907. (10) 
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From the above results it may be safely concluded that the first eigenvalue of Eq. (1) 
is \, = 723.907, correct to at least six figures. 

It is theoretically possible to obtain not only the first, but also the higher, eigen- 
values by this method. However, meaningful results can only be obtained for the higher 
ones if a very high degree of accuracy in the coefficients exists. With the present figures, 
the second highest roots of Eqs. (9) and (10) are respectively, .276 and .313. The latter 
value is probably correct to at least two places; the corresponding value of a3 is 32.0. 

It may be added that this same equation has been treated by Purday [3], also using 
a series in the independent variable. The results obtained there are 

a; = 2.83, as = 32 


which agree with the values of the present analysis. 


REFERENCES 
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A USEFUL INTEGRAL FORMULA FOR THE INITIAL REDUCTION OF 
THE TRANSPORT EQUATION* 


By RICHARD L. LIBOFF (Army Chemical Center) 


A theorem pertaining to spherical harmonics** states that if an angle 8 is determined 
by two directions, whose coordinates on the surface of a unit sphere are (6, a) and 
(6’, a’), and F(8) is a function expandable in a series of spherical harmonics in argument 


(cos 8), 
, 21+ 1 
F(@) = D jE kiPi(cos 6), (1) 


then 
[ F(8)Y,(0, a) dw = k,Y,(0’, «’). (2) 


Y (6, w) may either be spherical surface harmonic, a tesseral harmonic, or a spherical 


harmonic. 
The theorem has been noted because of its applicability in diffusion theory, for 


problems which militate use of the transport equation. 
The transport equation for various types of scattering (e.g., neutron, gamma, light) 


may in general, be written, 


1 ; re | ee 
| 8: VNG, 8) = —N + 7, | PONE, s’) dw’ + S. (3) 


*Received April 13, 1955. The topics of this paper appear in greater detail, in a report soon to be 
published by these Laboratories. (CRLR #479). 
**Hobson, Spherical and Ellipsoidal Harmonics, Cambridge, 1931, p. 146. 
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Equation (3) states that the spatial rate of change of the distribution function 
N(r, s) at the position r, in the direction of the unit vector s, is due to three effects which 
appear on the right side of the equation. They are, in order of their appearance, first, 
radiation lost through absorption, and scattering out of the direction s; secondly radia- 
tion which is scattered from all 42 solid angle into the direction s. p(@) is the angular 
differential cross section for the specific process involved. Thirdly there is the source 
contribution S. It is assumed in the above equation that a scatter is not accompanied 
by a change in wave length. u is the total narrow beam absorption coefficient. 

When the solution of (3) is through expansion processes, Eq. (2) becomes a useful 
integral formula for reducing the integral of (3) to its equivalent summation. 

For instance, consider the problem which permits the unit vector s to be replaced 
by the ordinate angle @ (e.g., spherical symmetry and infinite plane source). In this case 
the distribution function may be expanded as 


2l + 1 


re a,(r)P,(cos 6’). (4) 





N¢r, 8’) = > 


Substitution of this expansion into the kernel of (3) will transform the original 
integral into a series of integrals, each term of which is similar to expression (2). From 
this we may write for the integral of (3), 


i+ 
» 4r 


k, is the Legendre coefficient of the expansion of p(@). 

Following this with the replacement of N(r, s) by its equivalent expansion (4) into 
(3), will reduce the original integro-differential transport equation to a system of differ- 
ential equations involving the sequence {a,(r)}, knowledge of which completely deter- 
mines N(r, s). 

For more complicated geometries, where expansions in tesseral, or spherical surface 
harmonics are called for, the formula may be used in like manner, reducing the integral 
term of the transport equation to its corresponding summation, whence usually, a 
reduced system of equations is easily derivable. 


i) 


— 


a,(r)k.P,(cos 6), (5) 


TWO REMARKS ON HEISENBERG’S THEORY OF ISOTROPIC TURBULENCE* 
By WILLIAM H. REID! (Trinity College, Cambridge, England) 


1. Introduction. The exact dynamical equation for the rate of change with time of 
the energy spectrum function E(k) in isotropic turbulence may be written in the form 
dE(k)/at = T(k) — 2vk°E(k), (1) 


where 7'(k) is the transfer function usually denoted by this symbol. The incompleteness 
of Eq. (1) is well known and, in the past, several so-called “physical transfer theories” 
have been proposed in which a further relationship between 7(k) and E(k) is postulated; 





*Received May 9, 1955. 
1Present address: Exterior Ballistics Laboratory, Aberdeen Proving Ground, Md. 
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the specific form of the relationship then follows from the particular mechanism of 
energy-transfer considered and from certain general dimensional considerations. 

The theory which has attracted the greatest attention is the one due to Heisenberg [3] 
and is based on the coneept of an eddy viscosity; in this theory 7'(k) and E(k) are related 
by an equation which may be written in the form 

T(k) = ~ 98 [E(k’)/k?]'?? dk’ | k’?E(k’’) dk’. (2) 
dk Jk /0 
In this equation, « is, by hypothesis, a constant of order unity. A number of writers 
have attempted to derive the value of « under widely differing conditions and, in par- 
ticular, existing results show a large variation of x with the Reynolds number of the 
turbulence. 

One of the purposes of the present note, therefore, is to show that these results are 
partially in error and that in fact x, or more precisely the ratio S/x, where S is the 
skewness factor of —du,/dz, is practically independent of the Reynolds number and 
thus to remove one cause for criticism of Heisenberg’s theory. 

2. The behavior of S/x for small values of the Reynolds number. Since S is directly 
related to the second moment of the transfer function, one must first derive the explicit 
solution for 7(k). For small values of the Reynolds number, this can easily be done by 
expressing the solution as a power series in the Reynolds number. Thus, by assuming 


the usual type of initial period similarity 


E(k, t) = x v°7t"'?F(z) «and Tk, t) = «x ¥*?t°”?U(a), (3) 
where x = (vk*t)'””, Eqs. (1) and (2) become 
aF’(x) + (427 — 1)F(x) = 2U(z) (4) 
and 
U(2) = -2£ PF @)/2" ‘dz! | x! F(a") de’’. (5) 


From Eq. (5) it is clear that for small values of the Reynolds number U(z) is of higher 
order that F(x) and since R, , the Reynolds number of the turbulence usually denoted 
by this symbol, always occurs in the combination «xR, , we may write 


F(z) = >> F,(z)(kR,)” and U(x) = >> U,(x)(xR,)”. (6) 

n=2 n=3 
By substituting these series into Eqs. (4) and (5) and equating like powers of «R, , 
one obtains a sequence of equations for the determination of the functions F,(zx) and 
U,(zx). In particular, the first approximation to each of these functions can be readily 


obtained in the form 
F.(xz) = 2xe™**" (7) 


and 
U2) = 3)" 4 Bi (20 — (1 + 22’)e"**"] (8 
3\ 2) a 8 5 dz - \ \ ot Je ’ ) 


where —Ei(—z) is the exponential integral function defined by 


hin’ « [ dt. 


=z 
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These functions are shown graphically in Fig. 1. In passing, it may be noted that F,(z), 
unlike U;(x), is independent of the assumed form of the transfer theory. 








i | 
! 

2 

0.1- 

ok 1 l 
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Fic. 1. Curve 1: The ‘zero Reynolds number’ energy spectrum, F.(x); independent of transfer theory. 
Curve 2: The first approximation to the transfer function, 10U;(x); Heisenberg’s transfer 


theory. 


For small values of the Reynolds number, the ratio S/« is then given by 


Cae | Sth te + ORS. (9) 


K ‘ 0 


For the transfer function (8), the integral occurring in this equation has the value 


OY sia 3(15)'” (3 3) 
[ 2°Us(2) dx = “S— (5 — log 3 


and hence 
S,’x — 0.78 as R, — 0. (10) 
Experimentally, we know that S = 0.48 at the smallest Reynolds number for which 
experiments have been made [2, Fig. 6.3] and this makes x = 0.62. The numerical 
coefficient in Eq. (10) differs considerably, however, from the one given by Proudman [6]. 
Presumably, his result was obtained by expanding U;(x) in powers of x; in view of the 
somewhat peculiar behaviour of this function, the present method of using the exact 
form for U;(x) in closed form is obviously more reliable. 

3. The behavior of S/x for large values of the Reynolds number. For sufficiently 
large Reynolds numbers for which a universal equilibrium exists, the relation 


[ kE(k) dk 


iF =i 
| [ k? E(k) at 


S = 3 (39): ; (11) 
‘ 
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is exact; under these same conditions, one may therefore use Bass’ equilibrium spectrum 
to evaluate the behaviour of the two integrals occurring in this equation. This is the 
procedure used by Lee [5] who found that, as R, approaches infinity, 


S/« — 1.52. (12) 


In the course of the present investigation, however, it was found that when R, equals 
infinity, the behaviour of S/« is singular with a value which differs from the one given 
by Eq. (12) and this fact, in itself, may not be without interest. 

Thus, when R, equals infinity, the ratio S/x« can be determined by comparing the 
asymptotic form of the spectrum given by Kolmogoroff’s universal equilibrium theory 
with the corresponding form given by Heisenberg’s theory. For sufficiently small values 
of r, Kolmogoroff’s prediction for the double correlation function (see, for example, 
[1, p. 85}) 


2(u*)[1 — f(r)] = (4/58)**(er)””, (13) 


where (u’) is the mean square value of one component of the velocity, ¢ is the rate of 
viscous dissipation and f(r) is the double velocity correlation coefficient usually denoted 
by this symbol, leads to the equilibrium spectrum 


EQ) = 2 (ia) — (14) 
) = 81 (1/3)! \5S/ § , 


and this is to be compared with Heisenberg’s form of the equilibrium spectrum 
E(k) = (8/9«)?7e°k*”. (15) 


° a ome 2/3 5/3 = ° ° 
By equating the coefficient of e*k in these two expressions, one obtains 


S ] WP ee 2 ' 


and this then is the value of S/x« when R, equals infinity. 

From a physical point of view, the correct limit is of course the one given by Eq. (12), 
and when this value is used in conjunction with the experimentally determined value 
of S for large values of the Reynolds number, about 0.30 [1, Fig. 6.3], we obtain the 
value x = 0.20. 

4. Conclusions. The implication of this discussion then is that the ratio S/x varies 
between 0.78 and 1.52 as the Reynolds number varies from zero to infinity and, from 
the nature of the theory, it is reasonable to suppose this variation to be monotonic. 
Furthermore, when the experimentally observed variation of S with the Reynolds 
number is taken into account, it is found that « then varies between about 0.62 and 
0.20, which neatly bracket the value (0.45 + 0.05) suggested by Proudman [6] from his 
study of the correlation functions. It thus appears that the variation of « with the 
Reynolds number is relatively small and cannot be said to constitute a serious criticism 


of Heisenberg’s transfer theory. 
Acknowledgement. I am indebted to Dr. G. K. Batchelor for his helpful discussion 


of the present work. 
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NOTE ON LINEAR PROGRAMMING* 
By CARL E. PEARSON (Harvard University) 


Statement of theorem. Consider a set of m linear equations in n unknowns 7; ; 

Qq;Xji = ba ’ (1) 

where Greek indices range from 1 to m and Latin from 1 to n. The usual summation 

convention on repeated indices is used; a,; and b, are constants. Linear programming 

is a method of obtaining a solution (if it exists) of Eq. (1) satisfying in addition the 
requirements 

z; > 0, all i, (2) 

¢;x; = minimum, (3) 

where the c,; are constants. The fundamental theorem used in the simplex method of 

Dantzig (1) is that if one solution exists, then an equivalent solution can be found in 

which not more than m of the z; are non-zero; further, those columns of the matrix 

(a.;) which correspond to such non-zero (z,) will be linearly independent. The usual 

proof of this theorem (e.g. Ref. (2)) involves tedious geometrical considerations in 

n-dimensional space; it therefore seems worth-while to point out that a simple direct 


proof exists. 

Proof of theorem. Suppose (z/) satisfies conditions (1), (2), (3). Some—perhaps 
all—of these (x) will be non-zero; say for example that x; , 23 , 73 , are alone not zero. 
If firstly the corresponding columns (a,2 , Gas , @as) were linearly dependent, then a 
set of three constants K, , K; , Kz (not all zero) would exist such that 


A(Ke@a2 + Ks@a3 + Kedas) = 0, alla (4) 
for any arbitrary constant A. Then because of Eq. (4), the new set (z/’) defined by 


zi’ =2'— AK, fori = 2,3,6, 
(5) 


zi’ = 0 for other 7 


satisfies Eq. (1) and, for sufficiently small A, also Eq. (2). Clearly however c,zj’ < ¢;z} 
for appropriate A, unless 


C.K, + c3K3 + 6K, = 0. (6) 
*Received June 4, 1955. 








206 NOTES [Vol. XIV, No. 2 


Consequently, either the hypothesis of linear dependency has led to a contradiction 
so that the three columns are linearly independent, or alternatively Eq. (6) must hold. 
But if Eq. (6) holds, then 


C;2}' = C,4; for all A, 


1%; 
and A can be chosen so that at least one of x3’, 73’, 24’ vanishes, so that an equivalent 
solution involving fewer columns has been obtained; the same analysis may now be 


applied to the new solution. 
It follows that, eventually, an equally good solution utilizing only linearly independent 
columns will be obtained. Since not more than m columns can be linearly independent, 


the theorem is proved 
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HEAT FLOW IN A HALF SPACE* 


By WALTER P. REID{ (U. S. Naval Ordnance Test Station, China Lake, California) 


This heat flow problem has a more difficult boundary condition than usual. A formal 
solution is obtained by the combined use of Laplace and Fourier sine transforms. 

Let the region x > O have an initial temperature which is a function of x only, f(x). 
Assume that there is radiation of heat from the surface x = 0 to a finite slab, and from 
the other surface of the slab to surroundings whose temperature is a prescribed function 
of the time, g(t). The slab is assumed to be thin and of high thermal conductivity so 
that its temperature may be considered to be uniform throughout at any time. The 
heat exchange by radiation between adjacent surfaces will be taken to be proportional 
to the difference in temperatures of the two surfaces. Mathematically, then, the problem 


ill 


may be stated as follows: 
Ou(z, t) 0 ular, t > ra’ 
ee oe oi lor o> Gi >, (1) 
ot Ox 
u(x, 0) = f(x), (2) 
Au(O, t) 4 : 
on. tS a{u(O, t) — v(d)], t>OoO (3) 
Ox 
7 a1 { , 1 , 
h “ = a[u(0, t) — v(t)] — b(d) — gd], i> (4) 
Cl 
v(0) = V. (5) 


*Received June 13, 1955. 
tNow at Los Alamos Scientific Laboratory, Los Alamos, New Mexico. 
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In these equations u(x, t) is the temperature at some point in the region z > 0, 
while v(t) is the temperature of the slab. The quantities a, b, h and « are assumed to be 
positive constants. 

In the following both the Laplace transform and the Fourier sine transform will 
be used. This brings up the question of notation. A bar above the letter is sometimes 
used to represent a transform, or a change from capital letter to lower case, or a change 
in letter altogether. These methods are all satisfactory if just a single transform is used, 
but have disadvantages in the present case. In Eq. (16), for example, notation for three 
different transforms is needed. In this paper a transform will be indicated simply by 
changing the letter used in the argument of the function, as follows: 


ux, s) = [ e ‘u(x, t) dt, (6) 
ut, = / u(x, t) sin Ex dz. (7) 
The Laplace transforms of Eqs. (1), (3) and (4) are: 
2 
sul, s) — fz) = «ME (8) 
du(O, 8) _ _ 
Ox = alu(0, 8) v(s)], (9) 
hso(s) — hV = a[u(O, s) — v(s)] — bie(s) — g(s)]. (10) 


Hence 
fab — (a + b)0/dx + ahxd?/dx*® — hxd*/dx*Ju(0, 8) 
= abg(s) — ahf(0) + hf(0) + ahV. (11) 





Let 
w(x, t) = [ab — (a + b)d/dx + ahxd’/dx” — hxd*/dx*Ju(z, t). (12) 
Then 
dw(z, t) _—_-d*w(z, t) 
is” a (13) 
owls, 9) — w(x, = «Ue 9 | (14) 
sw(t, s) — w(k, 0) = xéw(0, 8) — Exwl, 8), (15) 
_ wk, 0) + xéw(0, 8) | 
w(é, 8) as s+ Ke (16) 


From Eq. (11): 
w(0, 8) = abg(s) — ahf(0) + hf’(0) + ahV. (17) 
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Take the Fourier sine transform of Eq. (12) and then eliminate the derivatives of u by 
integrating by parts. This gives 


A) 


w(t, t) = | u(x, t)[a(b — hxé*) sin Ex + E(a + b — Axé’) cos Ex] dx 


+ ahxtu(O, t) — hag O20 . (18) 
Let 
¢(~) = | f(x)[a(b — hxé’) sin tx + (a + b — het’) cos Ex] dz. (19) 
Then 
wt, 0) = P(&) + ahxtf(0) — hxtf’(0). (20) 
So 
fe — 00) + xéabyls) + xeahV , 
w(t, 8) = oe a (21) 
Hence 
w(t, t) = [o(t) + xtahV] exp (— xt) + xtab [ g(r) exp [— «é’(t — 7)] dr. (22) 
But 
9 
(23) 


w(x, i) == | w(t, t) sin Ea dé. 
Wd 0 
Therefore, from Eqs. (12) and (23), 
2f a(b — hxé’) sin Ex 


u(x, i) = - w(t, 5 > 
T | ? a’(b — hxé)° 








(a + b — hé)’ 


=: ee + b — ht) cos tx dé, (24) 


where w(E, ¢) is given by Eq. (22), and ¢() by Eq. (19). From Eq. (3) and (24) one may 


obtain v(?). 


NOTE ON THE SYMMETRICAL PROPERTY OF THE THERMAL CONDUCTIVITY 
TENSOR* 


By M. LESSEN (University of Pennsylvania) 


According to the Onsager reciprocal hypothesis, the thermal conductivity tensor 
among other similar properties of matter is symmetrical. In the following note it is 
demonstrated that for some cases of interest, the property of symmetry according to 
the Onsager hypothesis is not a necessary piece of information. 

Assuming a space filled with a continuous medium at temperature 7’, having a 








*Received June 15, 1955. 
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thermal conductivity of K,; such that the conducted flow of heat through an elemental 
area ds having associated direction cosines 1, is 


—(K;,,T lj ds. 


The net rate of heat flow into a volume bounded by a closed surface “s’’ then is 
[ (K,,T ,)l; ds. 
The surface integral may be transformed to a volume integral, 
[ Kut duds = [ Ku). do. 


It is therefore necessary to investigate the properties of the form 


(K,;T 3). ° 
Considering the case where 
K,; _ K;,(T), 
(K,;T ;).: ” K,;T 5; + aK pp, 


dT 


T;, and T ,7; are both symmetrical in 7 and j, therefore only the symmetrical portion 
of K,,; will matter in (K,,;7.;),; . 

Since the case of K;; = K,,;(T) applies to a large class of practical applications, it 
is important to note, that for this case the anti-symmetrical portion of K,; if it existed 
at all would not contribute to a first law of thermodynamic energy accounting. 


A CONVERSE TO THE VIRTUAL WORK THEOREM FOR DEFORMABLE SOLIDS* 


By W. S. DORN (Aircraft Gas Turbine Development Department, General Electric Company) 
AND 
A. SCHILD (Carnegie Institute of Technology) 


1. Introduction. Consider a continuous body occupying a volume V and bounded 
by a closed surface S.’ Any system of stresses o,; ,’ satisfying the equilibrium conditions 


for zero body forces 
o4;,;5 = 0, (1) 


Oi =~ OR, (2) 


*Received July 27, 1955. A. Schild’s participation was supported by a research grant from the 
National Science Foundation. 

1It is assumed that the body is simply connected and that the surface S is composed of a finite number 
of pieces of each possessing a continuously turning tangent plane. All of the functions will be assumed 
to possess as Many continuous derivatives in V and on S as are necessary for the theorems which will be 
used later. 

*The subscripts range over the values 1, 2, 3 and repeated subscripts will be summed over the entire 
range. Subscripts following a comma denote partial differentiation with respect to Cartesian coordinates 
2; ,C8-,0;,; = de; ;/ dx; . 
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everywhere in V, and a system of virtual displacements u,; in V, must together satisfy 
the equation of virtual work 


[ ose; dV = [ wsoun, dS, (3) 


where n; is the unit outward normal to the surface S, and where e,;; are the strains 


d 


derived from the displacements u; by 
é;, = #(u;.; + u;.,). (4) 


In this note the following converse to the theorem of virtual work is proved: 

If, for a symmetric tensor ¢;; given in V and for a vector u,; given on S, the virtual 
work equation (3) is satisfied for all equilibrium stresses o;; (i.e., for all o,; satisfying 
Eqs. (1) and (2)), then the ¢;; are compatible strains and are derivable, as in Eq. (4), 
from displacements u; which on S have the given boundary values. 

This theorem is similar to earlier results of R. V. Southwell’, and to some recent 
results of H. L. Langhaar and M. Stippes.* In contrast to the papers quoted, however, 
no stress-strain relations are assumed here and our theorem is not limited to linear 
elasticity theory. 


2. Stress functions. For each value of i = 1, 2, 3, Eq. (1) states that the divergence 


of a vector a; = o;; vanishes. Thus the vector may be expressed as the curl of another 
vector, so that 
6:15 = Aiinn » (5) 
Da. we mdhies « (6) 
Then Eq. (1) is identically satisfied. From Eq. (2) we have 
(Ain — Ajis).n = O. 
By the same argument as above, it follows that 
Aiin — Ajin = Biinmim » (7) 
Bijan = —Biinn = —Biinm - (8) 


Using the symmetries given by Eqs. (6) and (8), Eq. (7) may be solved for Ajj, , 
giving 
A jin _ 4(Baiim sa Binns + | e 


Therefore, by Eq. (5), 


0:5 = (Brim + Bimni),mn - (9) 
We now define 
Puiimn = (Basin + Bimni)- (10) 
Then 
re ne (11) 





’Proc. Roy. Soc. A 154, 4-21 (1936); Stephen Timoshenko 60th anniversary volume, The MacMillan 


Co., 1938, p. 211. 
‘J. Franklin Inst. 258, 371-382 (1954). 


1956] W. 8S. DORN AND A. SCHILD 211 


where, by Eqs. (8) and (10), P,;;, has the symmetries 
Pstten = —P isin = —Paini = Pi nni . (12) 


Thus any set of equilibrium stresses o,, are derivable, as in Eq. (11), from stress functions 
P iim Which satisfy the symmetry conditions (12). 

The six independent components of P,,;, can be expressed in terms of a symmetric 
second order tensor® 7',, which is a dual tensor of P,;im : 


Rw = des n€jeml niin ’ Pratt = eetinanT ws , (13) 


where ¢;,, is the usual completely skew-symmetric pseudotensor (¢,.3 = 1). The six 
components of P,,;,;,, or 7,, can also be identified with the stress functions x; , x2, Xs 
of Maxwell and y, , ¥2 , ¥; of Morera as follows:° 


X= —Pou323 =T1, Vi = 2Pai12 = —272; , 
Xe = —Pys:1 = T2: , Vo = 2Pio23 = —27 3: , (14) 
Xs = —Pisis = Ts ; Vs = 2P233, = —27 12 . 

For plane stress, the only non-vanishing stress function, x3 = —P1212 = 733 , reduces 


to Airy’s stress function. 
3. Compatibility equations. We shall now prove the first part of the theorem 
stated at the end of the Introduction, i.e., that the given stresses e;; are compatible. 
Expressing the equilibrium stresses o;; in terms of the stress functions P,,;, , Eq. (3) 
becomes 


| Pastmsmat dV = |e Pasim na dS. (15) 
Applying Green’s theorem twice to the volume integral, we obtain 
[ } er ee dV —= | SS A + Fetstltncus read P;; iam. m€in}; dS. (16) 


This equation must be valid for any equilibrium stresses, and thus for an arbitrary 
choice of the stress functions P,,;, . Let P.:;m vanish identically outside of a small 
volume surrounding an interior point P of V, and let P,;;, be essentially constant 
inside the small volume. Then the surface integral on the right side of Eq. (16) vanishes. 
Since the P,,,;,, at P are still arbitrary except for the symmetry conditions (12), it follows 
that 


Eig,mn — €nj,im — €im,in + Enm,iG = 0 (17) 


at any point P in V. By continuity, these equations are also valid on S. 
The equations (17) are the compatibility equations for strains and are necessary 
and sufficient conditions for the existence of a set of displacements U; in V such that 


i= (Ux; + U;,;). (18) 





5C. Weber, Z. f. Ang. Math. und Mech. 28, 193-197 (1948). 
6A. E. H. Love, A treatise on the mathematical theory of elasticity, 4th ed., Dover Publications, 
1944, p. 88. 
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It remains to prove that these displacements may be so chosen, that they take on 
the given boundary values on S. This will be done in the next section. 

4. Converse to the virtual work theorem. From Eq. (18), the (direct) theorem of 
virtual work for any equilibrium stresses follows by Green’s theorem: 


| esses dV = / Ue, a8. (19) 
Comparing with Eq. (3), we obtain 
[@—U)T.dS =0, Te = osm. (20) 


This equation must hold for any choice of tractions 7; on S which are obtainable from 
a distribution of equilibrium stresses o,;; in V. It is well known that any distribution of 
surface tractions 7’; which is in static equilibrium can be so obtained. 

Let Q and Q’ be any two points on S. Choose 7, to be identically zero outside of 
two small areas dS° and dS® which respectively surround the points Q and Q’. Inside 
these areas, let the tractions be essentially constant and such that the forces 


FS = TS dS°, F?’ = T?' dS*’ 
are equal and opposite and act in the direction of the line joining QQ’: 
Fo = —F? = X(zr% — 22), (A * 0). (21) 
This system of tractions is in static equilibrium. Equation (20) now gives 
[u; — Ude — (us — Ude-](xz? — x9’) = 0. (22) 
This equation is equivalent to the geometric statement that the (infinitesimal) 
displacement u; — U; leaves unchanged the distance between Q and Q’. Since this 


result holds for any two points on the closed surface S, it follows that u; — U,; must be 
a rigid body displacement: 


u; = U; + 0,;7; +; , (23) 
@;; = —w;; = const., @®; = const. 


Equation (23) holds on the surface S where the u; are given. We now use Eq. (23) 
to define throughout the volume V a displacement u, . It is then clear that this displace- 
ment takes on the assigned boundary values on S. From Eqs. (18) and (23) we also have 


é;; = 4(u;,; + u;,,). (24) 


This completes the proof of the converse to the virtual work theorem. 

5. Conclusion. The converse to the virtual work theorem might be used in some 
problems on deformable solids, with given surface displacements u; , to find among the 
members of an N parameter family ¢,;(z, , 22 , %3 } @1 , @ , *** , @y) Of strain functions 
the one which is ‘‘least incompatible’. 

One could choose a set of M equilibrium stresses a}; , --- , o:; and determine the 
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parameters a from the minimum principle 


M 
Peis ida be il oje;,;dV — [ won, as| 
A-l 


Alternatively, the parameters a could be determined from a minimax principle, such as 


Min,, ..-0» { Max, | / (> Buoit es av - ful B.o.! as}, (26) 


where the parameters 8 must satisfy }>4, 64 = 1. 


2 


(25) 


A NOTE ON LAMINAR AXIALLY SYMMETRIC JETS* 
By MARK BERAN (Wellesley, Mass.) 


Summary. It is shown that there is no stream function of the form y = rf(@), that is 
compatible with the complete Navier-Stokes equations, which represents a jet issuing 
from a small circular hole in an axially symmetric cone. 

The asymptotic velocity field of a laminar viscous jet is generally accepted to have 
a stream function of the form y = rf(@), corresponding to self-similar flow (Schlichting 
[1], Squire [2], and Yatseev [3]). The authors referred to have based their discussion on 
the fact that this assumption of self-similarity is compatible with both the boundary 
layer equations, and with the full Navier-Stokes equations. 

The purpose of this note is to establish a serious shortcoming of such models. It is 
shown that there is no continuously differentiable velocity field associated with a stream 
function of the form y = rf(@), which satisfies the Navier-Stokes equations and also 
adheres to a conical wall 6 = a > 0. 

Specifically, if Y = rf(@), then the velocity components in the r and 6 directions 
are respectively [4] 








a 
= E sin | dé’ (1) 
uy = | = | (2) 
*~ Lrsin 6° * 
The Navier-Stokes equations are equivalent to [5] 
f? = 4v cos Of — 2vsin oc — 2(c, cos’ 6 + c, cos 6 + ¢,) (3) 


for suitable constants c, , C2 , Cc; . We shall show that there is no solution of (3) which 
(¢) makes u, and wu, continuous for r > 0, and (iz) satisfies u,(a) = ue(a) = 0, for 0 


<ac<r. 
To show this, we also consider the differentiated form of (3), which is 


—f df _.d | 1d 4 
“aa eE - = : 4 
sin 646 ~ J — 28 075] sin 6 ag] — 2% 008 9 +e) (4) 
*Received June 28, 1955; revised manuscript received October 5, 1955. The work on this paper was 
partly supported by Contract N5ori-07634 with the Office of Naval Research. 
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The boundary conditions implied by conditions (7) and (72) are: 
For Eq. (4): (a) up = 0, u, finite, when 0 = 0. 
(b) u, = up = 0, when 6 = a. 


For Eq. (5): (c) us = 0, u, finite, du,/d0 finite, when @ = 0.* 
These yield respectively the equations: 


tate =0, (5) 
c, cos a +c, cosat+e, = 0, (6) 
2c, +c, = 0. (7) 
These have c, = c, = ¢; = 0 as their only solution. 


As Squire [5] has shown, the general solution of Eq. (4) with c, = c. = cs; = Ois: 
2v sin’ @ 
es ere ®) 
at! cos 6 


where a is an arbitrary constant. 
Referring again to the boundary condition u,(a) = ue(a) = 0, we see that there is 


no finite value of a(a = © yields a satisfactory but trivial solution) which satisfies this 
boundary condition, no matter what value of a is chosen. 
Thus we have shown that there is no non-trivial solution of the form y = rf(@) 


that is compatible with the Navier-Stokes equations and the boundary conditions (7) 
and (72). 

The author wishes to thank Prof. Garrett Birkhoff for suggesting the problem and 
for his helpful advice. 

REFERENCES 
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Schlichting, Grenzschichtstheorie, Karlsruhe 
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[3] V. L. Yatseev, Zh. eksp. teor. fiz. 20, 1031-34 (1950) 
{4] S. Goldstein (ed.), Modern developments in fluid dynamics, vol. 1, Oxford, 1938, pp. 104 and 115 
[5] H. B. Squire, Quart. J. Mech. Appl. Math. 4, 321-29 (1950) 


*It is assumed that when these boundary conditions are substituted in Eqs. (4) and (5) that the 
limit @ — 0 is taken, since @ = 0 is a singular point in the spherical polar coordinate system. 


HEAT CONDUCTION IN SEMI-INFINITE SOLID IN CONTACT 
WITH LINEARLY INCREASING MASS OF FLUID* 


By C. C. CHAO ann J. H. WEINER (Columbia University) 


Introduction. Problems of transient heat conduction in which the surface of a 
solid is in contact with a well-stirred fluid have been the subject of numerous investi- 
gations.’ In all previous cases studied, the mass of the fluid has been considered constant. 
However, it is sometimes of interest to know the temperature in the solid and fluid 





*Received August 9, 1955. 
1A review of previous work is found in [1], pp. 16-17. 
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while the latter is being poured onto the solid. This is so, for example, in the analysis 
of the filling stage of a casting process. 

In this paper, the problem of a semi-infinite solid in perfect contact with a well- 
stirred fluid whose mass is initially zero and increases linearly with time is considered. 
The principal boundary condition involved has variable coefficients and the solution of 
the problem by the Laplace transform technique requires the solution of a differential 
equation with the Laplace transform parameter as independent variable. Numerical 
results are presented. 

Problem. The problem described above is formulated mathematically as follows. 
(Property values are taken independent of temperature; the initial temperature of the 


solid is taken as zero.) 


07u Ou 
Ka 3? 2 > G; t> 0, (1) 
limu(z,t)=0; 2«>Q0, (2) 
1-0 
lim u(x, t) = 0; t> 0, (3) 
, Ou 0 = i _ 
K ee me 2 (tu) — V | at z= 0, t>0, (4) 


where u(z, ¢) is the temperature in the solid at distance x from the surface at time ¢ 
after the start of pour, K is the thermal conductivity of the solid, « its thermal diffusivity, 
V is the pouring temperature of the liquid, c is the specific heat of the liquid, and m 
is its (constant) mass rate of flow. 

Equation (4) is based on the assumptions that the fluid is well-stirred and is in 
perfect contact with the solid, that is, that the entire mass of liquid is at the temperature 
of the solid surface. It is derived from a heat balance on the liquid. 

Solution. The solution of the problem is obtained by use of the Laplace transform. 

Let U(z, p) = Jo €°”'u(z, t) dt denote the Laplace transform of u(z, t). Then the 


Laplace transforms of Equations (1)-(4) are: 








277 r 
oe, 2o@ (5) 
Ox K 
lim U(x, p) = 0, (6) 
. aU au . V 
[> = — mel » 32 +7], z= 0. (7) 


The transform of Eq. (4), Eq. (7), is obtained by using the operational properties 
of the Laplace transform ((2], operations 3 and 8, p. 294). 
The solution of Eq. (5) which satisfies Eq. (6) is: 


U = A(p) exp [—(p/x)'z]. (8) 
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Substitution of Eq. (8) into Eq. (7) leads to the following differential equation for A(p): 
dA 8A V (9) 


where s = K/mcex’”’. 
The solution of Eq. (9) is 


V 217 1/2\-1 1/2 * 1/2 1/2 
A(p) = -* 43°V[(2sp'’")"" + exp (2sp'’")Ei(—2sp’")] + C exp (2sp’"), — (10) 


where Ei(—z) = {2 y~'e~” dy is the exponential integral. It is seen that the arbitrary 
constant, C, must be zero if A(p) is to have a sectionally continuous inverse transform 


of exponential order. 
Therefore, substitution of Eq. (10) with C = 0 in Eq. (8) yields the transform of 


u(z, t) as 
U(z, p) = Vip"' — 4s°[(2sp'*)-* + exp (2sp'*)Ei(—2sp'”*)]} exp [—(p/x)'*z]. (11) 

The following transform pair is known ((4], p. 268, pair 31): 

L~*{(2s)~' + p exp (2sp)Ei(—2sp)] = (t + 2s8)~?. 

Therefore, 
L~"{{(2s)"* + p exp (2sp)Ei(—2sp)] exp (—apx'””)} = (t — xx”? + 28)", t> xx”. 
= 0, ¢<2.°™. 
The following property of Laplace transforms is known ((1], p. 243): If 

L~"[W(p)] = w(d), 


then 
L"(W(p'”)p"'””] = (wt)? | w(y) exp (—y’*/4t) dy 


therefore, 
L~ [exp (2sp'””)Ei(—2sp'’”) + (2sp'’”)~'] exp [—(p/x)'”z] 
= (rt)? [| (y — 2a? + 28)? exp (—y°/41) dy, 


= (te) [fe — [2/2 ut)"7] + 7}? exp (—2*) dt. 


Use of the above transformation pair in Eq. (11) together with the known pair ((1], 
p. 381, pair 8) 


L~"{p7* exp [—(p/x)'z]} = erfe [x/2(xt)'”*] 
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in Eq. (11) yields 


u(x, t) = v| exte [x /2(«t)'’?] 


- (12) 
- ast | ME — [2/2(xt)] + 807}? exp (—€) at]. 
z/2(«t)?/2 
The above analysis has been purely formal. It may be verified, however, by direct 
substitution that the above function satisfies Eqs. (1)-(4) defining the boundary value 
problem. The interchange of differentiation and integration required in this verification 
may be readily justified. 
Numerical results. The temperature of the liquid, which is equal to u(0, #4), is of 
particular interest. By use of integration by parts, it may be put in the following form: 


u(0, t) = v4 — 2s(rt)-'”? + 28°t* — 48°(rt*)~'”? [ (¢ + st-’”*)~* exp (—#’) ax. (13) 


Numerical values for the latter integral have been tabulated by Goodwin and Staton 
[3]. However, in the evaluation of Eq. (13), it was found necessary to carry their asymp- 
totic expansion further in order to obtain sufficiently accurate results. A graph of 
u(0, t)/V as function of t/s’ is shown in Fig. 1. 
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ERRORS IN ASYMPTOTIC SOLUTIONS OF LINEAR 
ORDINARY DIFFERENTIAL EQUATIONS* 


By F. W. J. OLVER (National Physical Laboratory, Teddington, England) 
In a recent paper [1] R. L. Evans has described a general method for estimating 
errors in asymptotic expansions of solutions of ordinary linear differential equations, 


the basis of which is as follows. 
Let 


y(2) ~  ele)a~* (4) 


be the formal asymptotic expansion for large x of a solution y(x) of the differential 
equation 


Ly) = 2d play*” = 0, (3) 
vy=0 
in which 
p(z)= > b,,2% [v = 0,1, --- , m and each integer n(v) is finite]. 


The error in question is the difference between y(x) and the sum of the first N terms 
in its asymptotic expansion, and is given by 


v(x) = y(z) — u(z), (6) 
where 
ulz) = ¥cla)a’*. (5) 


It satisfies a linear differential equation of the form 
M(v) — Le) = Ly), (11) 


where MM is a linear differential operator of order n whose coefficients can pe found from 
the p,(z), and whose term in the nth derivative cancels in (11) with that of L(v). Ac- 
cordingly (11) is of order n — 1. In particular, if n = 2 equation (11) is of the first order 
and yields immediately an indefinite integral for v(x) from which bounds for v(x) can be 
obtained. 

The purpose of this note is to point out that Evans’ result is incorrect. No non-trivial 
differential equation of the form (11) can exist. If it did, the relation y(z) = v(x) + u(z), 
obtained from (6), would enable us to assert that the general solution of any nth order 
differential equation of the form (3) can always be made to depend on the solution of a 
non-homogeneous differential equation of order n — 1. In particular, we would always 
be able to express the general solution of any second-order equation of the form (3) as a 
finite combination of indefinite integrals of elementary functions. 





*Received August 16, 1955. 
{The notation and equation numbers are the same as in [1]. 
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This conclusion is borne out by an example given in [1] in which the method is applied 


to the equation 
2 


7 om 1 i 1 nm = 
y +( 2+)y (1+) =0, 


with solutions e*J,(x), e*K,(x). It leads to the result 


N-1 


y = >> c(a)x™** + (constant) x 2°*” é | Saat te 


a=0 


where the constant depends on N but not on x. This suggests that the derivatives of 
a**I,(x) and x«“’’K,(x) can be expressed as finite combinations of elementary functions, 


which is incorrect. 
If the analysis given in [1] is examined, it is found that the relations given between 
Eqs. (9) and (10) do not follow from (9) and (6). It would appear that the correct forms 


of (9) and these relations are given by 


¢o=p—-N, wz) ~ YO Cex, 


a=0 


in which event (7) transforms into itself with c, p replaced by C, o respectively; in 
particular the equation between (10) and (10a) is to be replaced by 


{A.(o, — a — 2) + A;}C.(a + 2) + A;C.(a) 
+ {(¢. —a — 1)(o, — a — 2) + Axo, — a — 1) + A, JC + 1) = O. (A) 


In constructing the differential equation for v(x) which “‘corresponds”’ to (A) it must be 


remembered that C(— 1), C(— 2), --- , C(— N) are non-zero, accordingly this equation 
is non-homogeneous. Its correct form is in fact given by 
Lv) = —Liy), 


a result which is otherwise immediately obtainable from (3) and (6). The proposed 


scheme is accordingly nugatory. 
This note is published with the permission of the Director of the National Physical 


Laboratory. 
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Tables of integral transforms. Based, in part, on notes left by Harry Bateman. Compiled 
by Staff of the Bateman Manuscript Project. McGraw-Hill Book Co., Inc., New 
York, Toronto, London, 1954. Volume I, xx + 391 pp., $7.50. Volume II, xvi + 451 


pp., $8.00. 


The two volumes in question arise from the Bateman Manuscript Project at the California Institute 
of Technology under the support of the Office of Naval Research. Prepared in part from notes left by the 
late Harry Bateman by A. Erdelyi, W. Magnus, F. Oberbettinger, and F. G. Tricomi, these two volumes 
contain transform tables arranged as follows: Volume J. Fourier Sine, Cosine and Exponential Trans- 
forms 118 pp., Laplace Transforms 173 pp., Meilin Transforms 59 pp. Volume IJ. Bessel Transforms. 
Hankel Transforms of the form se f(x)J, (zy) (zy)? dz, 83 pp., and the same for 1/»(zy), K,(z, y) and 
H,(z, y), 75 pp., Kontrovich-Lebedev Transforms, 3 pp. Miscellaneous Transforms. Fractional Integral, 
Stieltjes and Hilbert Transforms, 78 pp. Integrals of Higher Transcendental Functions. Orthogonal 
Polynomials, Gamma functions and related functions, Legendre functions, Bessel functions and Hyper- 


geometric functions, 152 pp. 
Tables are very well arranged. The notation is explained in appendix to each volume. The printing 


job is really excellent. 
Roun TRUELL 


‘ 


Proceedings of the International Conference of Theoretical Physics. Kyoto and Tokyo, 
September, 1953. The Organizing Committee, International Conference of Theoretical 
Physics, Science Council of Japan, Ueno Park, Tokyo, 1954. xxviii + 942 pp. $10.00. 
(postage $1.00 extra). 


This book of more than nine hundred pages is a detailed account of a conference held in Japan in 
September 1952 which covered almost all of the currently interesting topics of theoretical physics. No 
attempt will be made to mention all of the authors who contributed to the conference. In fact it is un- 
fortunately impossible even to list the titles of all of the papers. There were about 125 papers presented. 
Because all of the topics covered were represented by the recognized leaders in these fields the relative 
emphasis in various fields is indicated below by the number of papers presented and the number of pages 
taken up in the proceedings. The main topics were: I. Field Theory and Elementary Particles. Papers and 
discussion concerned with Field Theory, Cosmic Rays and V-particles, Pions, Coupling Theory, and 
Nuclear Forces occupy approximately one third of this book. Forty five papers were presented in this 
area of physics with the result that these proceedings contain probably the best summary available of 
recent work. II. Nuclear Physics. Nuclear Reactions. Shell Structure and Beta Decay were the subject 
matter of a dozen papers and eighty pages of these proceedings. III. Statistical Mechanics. The Theory of 
Polymers, Liquids, Transport Phenomena, Irreversible Processes and General Methods of Statistical 
Mechanics were topics covered in twenty one papers and 180 pages of the proceedings. 

IV. Molecules and Solids. Dislocations, Molecular Theory, Metals, Electron Theory of Intrinsic 
Magnetization, Antiferro- and Ferrimagnetism, Magnetic Resonance, Dielectrics, and color centers were 
the subjects in that section of the conference dealing with solid state theory. Thirty seven papers cover 
311 pages of the proceedings. 

V. Liquid Helium and Superconductivity. Theories of Liquid Helium and Superconductivity are dis- 
cussed in ten papers covering about sixty pages of the proceedings.—This conference and the proceedings 
appear to have made a large contribution by putting into print an outline of the present state of theo- 


retical physics as it is seen by a large number of physicists from all countries. 
Roun TRUELL 
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Introduction to theoretical mechanics. By Robert A. Becker. McGraw-Hill Book Co., 
Inc., New York, Toronto, London, 1954. xiii + 420 pp. $8.00. 


This is a text book on a level appropriate for junior or senior year courses. It is based on a course for 
students in engineering physics, and the choice of topics, examples, and problems has a flavor of physics 
rather than of engineering. Some 80 worked and 400 unworked problems are given. Vectors are used 
throughout, but tensor or matrix concepts are not mentioned. The chapters on plane and general motion 
of a rigid body seem especially well written and complete, considering the elementary level of mathe- 
matics used. The chapter on vibrating systems, however, seems inadequate since it is confined mainly 
to properties of normal coordinates, and the standard techniques of combining normal modes are not 
introduced. The principle of virtual displacements appears not to be given the importance it deserves as 
an analytical tool in both statics and dynamics (in deriving Lagrange’s equations, for example). Apart 
from these minor criticisms the reviewer was well impressed by this as a carefully written and reliable 


basic text. 
P. S. Symonps 


Geometrical mechanics and de Broglie waves. By J. L. Synge. Cambridge University 
Press, 1954. vi + 167 pp. $4.75. 


This monograph presents a rather unique amalgamation of several disciplines from the author’s 
unusual collection of interests. The result of this synthesis is an approach to relativistic particle dynamics 
based upon Hamilton’s methods of geometrical optics, the Minkowski geometry of special relativity 
and the particle-wave concept of de Broglie. 

After a short introductory chapter, the author devotes about half of the book to a generalization of 
the Hamilton methods to the four-dimensional Minkowski space-time. Although the methods are 
described as Hamilton’s methods of optics, it is pointed out that these differ only formally from Hamilton’s 
analysis of classical dynamics. The former, however, is in a form more easily adapted to the generalization 
in question. In this generalization, the rays of optics (trajectories of dynamics) become time-like world 
lines of particles. The waves of optics (2-dimensional surfaces) become 3-waves in the Minkowski space 
describing the history of a 2-wave. 

One of the very interesting results of this procedure is that two velocities, a ray velocity and a wave 
velocity emerge from the theory in a natural way. These are interpreted to correspond respectively to 
the particle and wave velocities of de Broglie waves. These two velocities are obtained, however, solely 
from geometrical considerations without any discussion of wave interference or other concepts cus- 
tomarily involved in the association of particle velocity with a group velocity. 

Not until chapter IV does the author attempt to assign a phase to these 3-waves. He does so then 
only for the purpose of incorporating into the theory a “primitive quantization” in much the way that 
one obtains physical optics from geometrical optics by the assignment of a phase to the wave front. 
Several special applications of this relativistic quantum mechanics are considered, including the hydrogen 
atom quantization, the Zeeman effect and interference of a particle from two holes. 

The final chapter is concerned with the further generalization of the methods to the two body 
problem in an 8-dimensional product of two Minkowski spaces. 

The author points out that his primary concern is to develop a “coherent mathematical theory”’; 
that physical interpretation is considered of secondary importance. Indeed the theory is not all inclusive 
from a physicist’s point of view, for there is no obvious way of incorporating into the theory the concepts 
of spin, or “second quantization” and even the “first quantization” according to this scheme has been 
tested on only a few problems. Where the theory is known to give correct or good approximate results, 
however, this approach to the physical problem furnishes a very elegant geometrical interpretation. 

The purpose of the theory is well stated, the postulates are well defined and the development proceeds 
in a systematic fashion. Although this reviewer thinks that the book gives an excellent presentation of 
the ideas that it attempts to convey, he fears that it will not receive the circulation it deserves because 
it does not seem to belong specifically to any of the many narrow channels of current scientific interest. 
The physical principles underlying the theory have been known since 1928. If this book could have 
been published then, it would undoubtedly have been an immediate sensation. 

Gorpon F. NEWELL 
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Modern developments in fluid dynamics. Two Volumes. Edited by L. Howarth. Oxford 
University Press, London, 1953. xvi + 875 pp. $17.00 per set. 


These two volumes were developed as companion volumes to ‘‘Modern developments in fluid dynamics” 
edited by S. Goldstein, which appeared in 1938. These latter books have become indispensable to all 
those working in fluid dynamics, since they bring together in a concise and accurate form the accumu- 
lated understanding of fluid mechanics phenomena as of the time of their publication. The present 
companion volumes on high speed flow are also sponsored by the fluid motion panel and attempt to 
similarly summarize developments in this latter field. 

The publication of two new volumes instead of a revision of the older volumes is possible because 
much of the new work in fluid mechanics which grew out of wartime problems has been in the field 
of high velocities. However, any publication in this latter field must be carried out under the disadvantage 
of attempting to describe the important phenomena without violation of security restrictions which 
arose during the war and which have only been partially relaxed. As in the previous volumes, a group 
of distinguished authors have contributed sections to the new work. 

Roughly, the first volume is devoted to discussions of theory and the second volume to questions 
of experimental results. The first two chapters are devoted to a general introduction to the various high 
speed phenomena and to the equations in both vector and tensor notation required for their descrip- 
tion. In these chapters, as throughout the book, the editor has attempted to make clear the distinction 
between flows in which the entropy is constant for each particle, so that its change of state is isentropic, 
and flows in which the entropy is constant everywhere. These latter are described as homentropic. 
While this distinction is clearly a significant one, the reviewer was sometimes confused by the words 
used and was left with the feeling that the change was somewhat of a fetish. 

Chapter 3 on the characteristics method is followed by Chapter 4 in which shock waves and their 
formation and various interactions are described. After thus treating fundamentals in some detail, 
Chapter 5 presents the few available exact solutions, vortex flow, radial flow, spiral flow and Prandtl- 
Meyer flow including their various applications. 

Chapter 6 and 7 treat the one-dimensional approximation with applications to nozzles and jets, and 
the various techniques available by use of the hodograph plane. Then comes several chapters dealing 
with approximate methods. The various linearization procedures and numerous of their airfoil, planform, 
tube, and jet applications are followed by available solutions to unsteady problems, primarily unsteady 
airfoil characteristics 

Chapter 10, the last one in the first volume, treats the general theory of boundary layers including 
skin friction, aerodynamic heating, heat transfer, stability, and boundary layer-shock wave interaction. 
The first volume contains 475 pages. 

The second volume begins with Chapter 11 devoted to experimental methods which is broken down 
into four sections. The first, wind tunnels and moving bodies, discusses general tunnel characteristics 
and supersonic nozzle design together with very brief statements about intermittent tunnels, shock 
tubes, and moving body tests. It appears to the reviewer that the latter two methods have been relatively 
neglected since they receive less than 2 pages of attention. 

The second section devoted to uncertainties and corrections arising in tunnel tests is followed by a 
section on measurements and their reliability, and finally, a section on visualization in which each of the 
common optical systems is treated in some detail. In Chapter 12, the experimental results available on 
flow past airfoils and cylinders are reported together with some comparisons with theoretical predictions. 
In this chapter, the lift, drag moment and their various derivatives are discussed with respect to variation 
with Mach number, Reynolds number, angle of attack, angle of sweep, section parameters and controls. 

Chapter 13, devoted to flow past bodies of revolution is devoted primarily to experimental and 
theoretical results on projectiles. The base pressure receives very considerable consideration. The range 
testing technique and the theory of the spinning shell are discussed. Chapters 11, 12, 13, which con- 
stitute three-quarters of the second volume are the chapters wherein restrictions on publication make a 
complete discussion of the problems difficult. The authors have attempted an important piece of writing 
using generally available information to bring out as clearly as possible the essential features of ex- 
perimental results in the high-speed aerodynamic field. 

The final Chapter 14 on heat transfer would better be considered a revision of what appears as 
the last chapter with the same title in the earlier volumes. Many sub-sections are identical to the previous 
volume, although many others have been rather completely rewritten, so as to bring them up to date. 
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In this heat transfer chapter, there is no attempt to limit the discussions to “high speed flow” although 
this phase of the subject has similarly been revised from the previous editions. 
These two volumes are, the reviewer believes, destined to take their position beside the earlier two 
volumes as the major publications in the field of fluid mechanics during the middle of our century. 
Howarp W. Emmons 


Vector and tensor analysis. By G. E. Hay. Dover Publications, Inc., New York, 1953. 

viii + 193 pp. $2.75. 

Many excellent books are available on both the vector and tensor ‘calculus and their applications- 
The present book contains no preface, and accordingly its purpose must be inferred from its contents. 
The book seeks to develop the essential parts of the vector and tensor calculus together with some dis- 
cussion of geometrical and physical applications. This, a formidable task to perform within a monograph, 
is made more difficult by reducing to a minimum the supposed background knowledge of the reader. 
The greater part (156 pages) of the book is devoted to the vector calculus, covering the usual material 
up to Green’s and Stokes’ theorems, and including applications in geometry and rigid body dynamics. 
This part is well written, and a clear account is given. The remainder (37 pages) of the book is devoted 
to the tensor calculus with very brief mention of its applications. Here the coverage (this includes weighted 
tensors and covariant differentiation) is achieved only through extreme condensation of analysis, and 
the average reader would find it difficult to work through the probiems. 

The book contains a few incorrect statements and misprints. A partial bibliography is given but it 


is surprising to find no reference to several important standard texts. 
H. G. Hopxins 


Proceedings of the Eastern Joint Computer Conference, held by the Joint Computer 
Committee of the Association for Computing Machinery, American Institute of 
Electrical Engineers, and the Institute of Radio Engineers, Inc., in Philadelphia, 
December 8-10, 1954. American Institute of Electrical Engineers, New York, 1955. 


92 pp. $3.00. 


The theme of this conference was the design and application of small digital computers. The Pro- 
ceedings Volume commences with an introductory talk by C. W. Adams, who defines small digital 
computers by means of examples. He identifies two classes: (1) Those built up from punch card equip- 
ment, which he calls “‘big-little computers’’, and (2) Those built up around a magnetic drum, which he 
calls “‘little-big computers’’. (The significance of these names for the two classes would perhaps be 
clearer if the hyphen in each case were moved forward by a word.) Computers of the second class typically 
have an internal storage capacity of between 10,000 and 40,000 decimal digits. 

The papers presented at such a conference fall into several recognizable categories. First, there is the 
category of survey papers, which discuss available equipment or current trends in development. In the 
present case, this category is represented by a good survey of currently available small digital computers 
presented by A. J. Perlis. Second, there are the papers devoted to engineering design problems. For 
various reasons, some of which are quite obvious, these have a tendency to be somewhat vague when 
presented by engineers working for private industry. The volume under review contains three or four 
papers of this type covering the following subjects: pulse packing density, magnetic core circuits, and 
automation of information retrieval. 

A third recognizable type of paper is the company report which describes a specific proprietary piece 
of equipment. As is usual at computer conferences held by the engineering societies, this category is well 
represented in the current volume. 

Finally, there are the contributions of those interested in the applications. In the present conference 
the emphasis here was mainly on business problems. However, there are two papers in the Proceedings 
Volume with some mathematical content. The first is a paper by H. M. Gurk and Morris Rubinoff on the 
numerical solution of differential equations. This paper contains some interesting stability charts for the 
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finite difference analogue of the equation ¢ = \z where A is complex. The second mathematical paper is a 
simplified discussion of the use of digital computers in optical ray tracing, by N. A. Finkelstein. 

The reviewer found this volume to be a relatively interesting one of its kind and the editorial work 
was excellent. But this fact did not deflect the reviewer from the pursuit of a project which he has had 
in mind for some time. The project is to found a Society for the Prevention of the Publication of Proceed- 
ings of Symposia. The holding of symposia is a fine thing indeed, but the morning after (which sometimes 
takes place two years later) is a different matter altogether. The papers with substantial content pre- 
sented at a symposium should be, and almost always are, published elsewhere than in the Proceedings, 
with fuller detail. The papers which do not have substantial content should not be published at all. 

Membership in the new Society is free. Any member who can prove that he has been instrumental 


in the suppression of a symposium proceedings volume will be made a Fellow. 
J. H. Curtiss 


Table of Everett’s interpolation coefficients. By E. W. Dijkstra and A. Van Wijngaarden. 
Excelsior’s Photo-Offset, The Hague, 1955. $1.58. 


This table was constructed on the electronic computer ARRA of the Mathematisch Centrum in 
Amsterdam on request of Empresa Nacional Bazan de Construcciones Navales Militares in Madrid. 
It gives to seven decimal places the values of the coefficients of the second, fourth, and sixth central 
differences in Everett’s interpolation formula. 

This table is especially convenient for three reasons; a) Entries for the argument p are given at 
intervals of 0.0001 so that interpolation is rarely needed in the table itself. b) All the coefficients needed 
for any given interpolation are printed in the single line of the table, i.e., the coefficients for p and for 
1 — p are both given on the same line. c) The range extends from p = 0 to p = 1, so that it is not 
necessary to change the procedure for p > 1/2. As a result of b) and c) all entries appear twice in the 
table, once at p and once at 1 — p. 

The use of sixth order differences in Everett’s formula is equivalent to using eight equally spaced 
values of the function for the interpolation, or in other words, using a polynomial approximation of 


seventh degree. 
W. E. MILNe 


Tables of the cumulative binomial probability distribution. By The Staff of the Computation 
Laboratory. Harvard University Press, Cambridge, 1955. 1xi + 503 pp. $8.00. 


This, the first set of tables to be computed on the new Harvard Mark IV Calculator, gives values 
of the cumulative binomial probability distribution 
E(n, r, p) =>.2., p?(1 — p)*-2 n!/z!(n — z)! 
to five decimal places for integer values of r and n. The unique feature of the tables is that n ranges in 
varying steps from 1 to 1000 as compared with a previous high of 150. Since only integer r and n are in- 
cluded, the tables are not quite as complete for small n as the old tables by Pearson for the incomplete 
beta function which included half integer values as well. The incomplete beta function differs from 


E(n, r, p) only in a trivial way. 
The introduction includes seventeen illustrative applications in addition to a discussion of the 


analytic properties of the function. 
G. F. NEwELL 




















